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Abstract 


In  this  paper,  we  address  the  inverse  problem  of  source  reconstruction  for  the  difficult  case 
of  multiple  sources  when  the  number  of  sources  is  unknown  a  priori .  The  problem  is  solved 
using  a  Bayesian  probabilistic  inferential  framework  in  which  Bayesian  probability  theory 
is  used  to  derive  the  posterior  probability  density  function  for  the  number  of  sources  and 
for  the  parameters  (e.g.,  location,  emission  rate,  release  duration)  that  characterize  each 
source.  A  mapping  (or,  source-receptor  relationship)  that  relates  a  multiple  source  distri¬ 
bution  to  the  concentration  data  measured  by  the  array  of  sensors  is  formulated  based  on  a 
forward-time  Lagrangian  stochastic  model.  A  computationally  efficient  methodology  for  de¬ 
termination  of  the  likelihood  function  for  the  problem,  based  on  an  adjoint  representation 
of  the  source-receptor  relationship  and  realized  in  terms  of  a  backward-time  Lagrangian 
stochastic  model,  is  described.  An  efficient  computational  algorithm  based  on  a  reversible 
jump  Markov  chain  Monte  Carlo  (MCMC)  method  is  formulated  and  implemented  to  draw 
samples  from  the  posterior  density  function  of  the  source  parameters.  This  methodology 
allows  the  MCMC  method  to  jump  between  the  hypothesis  spaces  corresponding  to  different 
numbers  of  sources  in  the  source  distribution  and,  thereby,  allows  a  sample  from  the  full 
joint  posterior  distribution  of  the  number  of  sources  and  source  parameters  to  be  obtained. 
The  proposed  methodology  for  source  reconstruction  is  tested  using  synthetic  concentration 
data  generated  for  cases  involving  two  and  three  unknown  sources. 
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Resume 


Dans  cet  article,  nous  abordons  le  probleme  inverse  de  reconstruction  d’une  source,  concer- 
nant  le  cas  difficile  de  sources  multiples,  quand  le  nombre  de  source  est  a  priori  inconnu. 
Le  probleme  est  resolu  au  moyen  d’un  concept  inferenciel  de  probabilite  bayesienne  avec 
lequel  on  utilise  la  theorie  de  probabilite  bayesienne  pour  deriver  la  densite  de  probabilite 
posterieure  pour  le  nombre  de  sources  et  les  parametres  (ex.  :  endroit,  taux  d’emission, 
duree  d’emission)  qui  caracterisent  chaque  source.  Un  mappage  (ou  une  relation  source- 
recepteur),  faisant  la  relation  entre  la  repartition  des  sources  multiples  et  les  donnees  de 
concentration  mesurees  par  le  reseau  de  capteurs,  est  formule,  base  sur  le  modele  stochas- 
tique  Langrangian  futuriste  (forward-time)  On  y  decrit  une  methodologie  computationnelle 
efRciente  visant  a  determiner  la  function  de  vraisemblance  pour  le  probleme,  basee  sur  une 
representation  adjointe  de  la  relation  source-recepteur  et  realisee  en  termes  d’un  modele 
stochastique  Langrangian  passeiste(  backward-time).  Un  algorithme  computationnel  effi¬ 
cient  base  sur  une  methode  Monte  Carlo-chaine  de  Markov  (MCMC)  de  saut  reversible  a 
ete  formule  et  implemente  pour  tirer  des  echantillons  provenant  de  la  densite  de  probabilite 
posterieure  des  parametres  de  la  source.  Cette  methodologie  permet  a  la  methode  MCMC 
de  sauter  entre  les  espaces  d’hypothese  correspondant  aux  nombres  de  sources  differentes 
dans  la  repartition  des  sources  et,  par  consequent,  permet  d’obtenir  un  echantillon  tire  de 
la  repartition  posterieure  des  sources  integralement  conjointes  du  nombre  de  sources  et  des 
parametres  de  ces  sources.  Cette  methodologie  proposee  pour  la  reconstruction  des  sources 
est  testee  au  moyen  de  donnees  de  concentrations  synthetiques  generees  pour  des  cas  ayant 
deux  ou  trois  sources  inconnues. 
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Executive  summary 


Introduction:  An  increasingly  capable  sensing  technology  for  concentration  measurements 
of  chemical,  biological  and  radiological  (CBR)  agents  released  into  the  atmosphere,  either 
accidentally  or  deliberately,  has  fostered  interest  in  exploiting  this  information  for  detection, 
identification  and  reconstruction  of  agent  sources  responsible  for  the  observed  concentra¬ 
tion.  A  critical  capability  gap  in  current  emergency  and  retrospective  management  efforts, 
directed  at  terrorist  incidents  involving  the  covert  release  of  a  CBR  agent  in  a  densely 
populated  urban  center,  is  the  determination  of  the  number  of  sources  and  for  each  of 
these  sources  the  location  and  amount  of  agent  released  (source  term  estimation  problem) , 
following  event  detection  by  an  array  of  independent  CBR  sensors.  In  order  to  address 
this  capability  gap,  an  additional  component  for  source  reconstruction  has  been  included 
in  the  Chemical  Biological  Radiological  and  Nuclear  (CBRN)  Research  and  Technology 
Initiative  (CRTI)  project  02-0093RD  (Advanced  Emergency  Response  System  for  CBRN 
Hazard  Prediction  and  Assessment  in  the  Urban  Environment).  This  component  has  has 
been  incorporated  in  CRTI  02-0093RD  under  the  auspices  of  the  Public  Security  Techni¬ 
cal  Program  (PSTP)  project  “United  States/Canada  Collaborative  Projects  on  Science  and 
Technology  in  Urban  Transport  Modeling  Related  to  Homeland  Security”  and  funded  under 
the  CRTI  Supplemental  Punding  Program  for  existing  projects,  with  Defence  R&D  Canada 
-  Suffield  as  the  principal  investigator. 


Results:  In  this  report,  we  address  the  problem  of  source  term  estimation  for  the  difficult 
case  of  multiple  sources  when  even  the  number  of  sources  is  unknown  a  priori.  A  new 
methodology  was  developed  to  solve  this  problem.  This  methodology  provides  simultaneous 
estimates  for  the  number  of  sources  and  for  the  parameters  (e.g.,  location,  emission  rate, 
release  duration)  that  characterize  each  source.  The  proposed  methodology  for  source 
reconstruction  has  been  tested  successfully  using  synthetic  concentration  data  generated 
for  cases  involving  two  and  three  unknown  agent  sources. 


Significance  and  Future  Plans:  The  algorithm  developed  in  this  report  may  be  used  to 
interpret  agent  concentration  measurements  obtained  from  an  array  of  sensors  in  order  to 
estimate  unknown  source  characteristics  and  quantify  the  uncertainties  in  this  estimation. 
The  technique  may  be  used  to  infer  the  location,  emission  rate,  and  duration  of  putative 
agent  release(s)  and  this  information  can  be  subsequently  used  to  predict  future  agent 
transport  in  the  atmospheric  environment.  In  summary,  the  algorithm  provides  a  method¬ 
ology  for  the  fusion  of  CBR  sensor  data  with  model  predictions  of  agent  dispersal  in  the 
atmosphere  to  provide  a  more  complete  situational  awareness  in  the  operational  environ¬ 
ment.  Future  plans  include  testing  the  proposed  Bayesian  inference  algorithm  for  multiple 
sources  using  real  concentration  measurements  obtained  from  a  cooperative  FUsing  Sensor 
Information  from  Observing  Networks  (FUSION)  Field  Trial  2007  (FFT  07)  that  will  be 
conducted  at  U.S.  Army  Dugway  Proving  Ground  (DPG)  in  September  2007.  This  field 
experiment  is  conducted  under  Technical  Panel  9  (TP-9)  of  The  Technical  Cooperation 
Program  (TTCP)  Chemical,  Biological  and  Radiological  Defense  (CBD)  Group. 


Yee,  E.  (2007).  Bayesian  Inversion  of  Concentration  Data  for  an  Unknown  Number  of 
Contaminant  Sources.  (DRDC  Suffield  TR  2007-085).  Defence  R&D  Canada  -  Suffield. 
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Sommaire 


Introduction:  Une  technologie  de  detection  de  plus  en  plus  efficace  a  mesurer  la  concen¬ 
tration  d’agents  chimiques,  biologiques  et  radiologiques  (CBR)  liberes  dans  Tatmosphere, 
soit  accidentellement  soit  deliberement,  a  suscite  I’interet  concernant  I’exploitation  de  cette 
information  pour  la  detection,  I’identification  et  la  reconstruction  de  sources  d’agents  ayant 
cause  la  concentration  observee.  Diminuer  I’ecart  crucial  entre  la  capacite  des  efforts  a 
repondre  aux  urgences  et  celle  de  leur  gestion  retrospective  dirigee  centre  des  incidents  ter- 
roristes  comprenant  la  dissemination  secrete  d’agent  CBR  dans  un  centre  urbain  densement 
peuple,  consiste  a  determiner  le  nombre  de  sources  et  pour  chacune  de  ces  sources,  I’endroit 
et  la  quantite  d’agent  dissemine  (probleme  d’estimation  en  termes  de  sources)  qui  suit  la 
detection  d’un  evenement  au  moyen  d’un  reseau  de  capteurs  CBR  independants.  Pour  etre 
en  mesure  de  diminuer  cet  ecart  entre  les  capacites,  une  composante  complement aire  de  la 
reconstruction  des  sources  a  ete  indue  dans  le  Projet  02-0093RD  de  I’Initiative  de  recherche 
et  de  technologie  chimique,  biologique,  radiologique  et  nucleaire  CBRN  (IRTC)  (Systeme 
d’intervention  d’urgence  avance  pour  la  prediction  et  revaluation  des  risques  CBRN  en  mi¬ 
lieu  urbain).  Cette  composante  a  ete  incorporee  dans  le  Projet  02-0093RD  de  I’lRTC  sous 
I’egide  du  projet  <C  Projets  communs  entre  le  Canada  et  les  Etats-Unis  sur  la  science  et 
la  technologie  de  la  modelisation  des  transports  en  commun  relies  a  la  Securite  interne  ^ 
du  Programme  technique  de  securite  publique  (PTST)  qui  est  finance  par  le  Programme 
de  financement  supplement  aire  de  I’lRTC,  pour  les  projets  existants,  avec  R&D  pour  la 
defense  Canada  -  Suffield  pour  instigateur  principal. 


Resultats:  Dans  cet  article,  nous  abordons  le  probleme  inverse  de  reconstruction  d’une 
source,  concernant  le  cas  difficile  de  sources  multiples,  meme  quand  le  nombre  de  source  est 
a  priori  inconnu.  On  a  mis  au  point  une  nouvelle  methodologie  pour  resoudre  ce  probleme. 
Cette  methodologie  procure  des  estimations  simultanees  pour  le  nombre  de  sources  et  leurs 
parametres  (ex.  :  endroit,  taux  d’emission,  duree  d’emission)  caracterisant  chaque  source. 
Les  essais  sur  cette  methodologie  proposee  ont  produit  des  resultats  positifs  pour  la  recon¬ 
struction  des  sources  au  moyen  de  donnees  de  concentrations  synthetiques  generees  pour 
des  cas  ayant  deux  ou  trois  sources  inconnues. 


La  Portee  des  Resultats  et  les  Plans  Futurs:  L’algorithme  mis  au  point  dans  ce  rap¬ 
port  pent  etre  utilise  pour  interpreter  les  mesures  de  concentration  d’agent  obtenues  d’un 
reseau  de  capteurs  en  vue  d’estimer  les  caracteristiques  d’une  source  inconnue  et  de  quanti¬ 
fier  les  incertitudes  contenues  dans  I’estimation.  La  technique  peut  etre  utilisee  pour  inferer 
I’endroit,  le  taux  d’emission,  la  duree  des  disseminations  de  I’agent  putatif  et  cette  infor¬ 
mation  peut  etre  consecutivement  utilisee  pour  predire  le  transport  futur  de  I’agent  dans  le 
milieu  atmospherique.  En  resume,  I’algorithme  procure  une  methodologie  pour  la  fusion  des 
donnees  des  capteurs  CBR  avec  des  predictions  d’un  modele  de  dispersion  d’un  agent  dans 
I’atmosphere  qui  visent  a  procurer  une  connaissance  plus  complete  de  la  situation  du  milieu 
operationnel.  Les  plans  futurs  consistent  a  evaluer  I’algorithme  d’inference  bayesienne  pour 
des  sources  multiples  utilisant  des  mesures  reelles  de  concentration  obtenues  en  cooperation 
de  I’information  fusionnee  des  capteurs  du  reseau  meteorologique  (ELISION)  des  Essais  sur 
le  terrain  2007  (EET  07)  qui  seront  conduits  sur  le  Polygone  d’essai  Dugway  de  I’arme 
americaine  en  septembre  2007.  Cette  experience  sur  le  terrain  est  conduite  sous  la  direction 
du  Panel  technique  9  (PT9)  du  Programme  de  cooperation  technique  (PCT)  du  groupe  de 
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Introduction 


The  ever  increasing  power  and  sophistication  of  today’s  chemical,  biological,  and  radio¬ 
logical  (CBR)  detection  technologies  are  opening  up  a  new  realm  for  the  detection  and 
monitoring  of  releases  of  these  toxic  agents  into  the  atmosphere.  Rapid  detection  and  early 
response  to  a  release  of  a  CBR  agent  could  dramatically  reduce  the  extent  of  human  expo¬ 
sure  (and  associated  morbidity  and  mortality  of  the  incident)  and  minimize  the  cost  of  the 
subsequent  clean  up.  In  the  event  of  a  CBR  incident,  the  assessment  of  the  damage  likely 
to  be  caused  by  the  release  is  a  problem  of  great  importance.  This  assessment  is  usually 
undertaken  using  a  predictive  model  for  the  mean  transport  and  turbulent  diffusion  of  the 
contaminant  through  the  atmosphere,  which  in  turn  provides  the  information  required  to 
determine  the  temporal  window  and  geographical  extent  of  the  hazard  zone  (or,  toxic  cor¬ 
ridor)  required  in  the  formulation  of  an  effective  response.  Unfortunately,  an  array  of  CBR 
sensors  by  itself  is  not  sufficient  for  this  task,  owing  to  the  fact  that  detection  of  a  toxic 
agent  plume  by  the  sensor  array  only  indicates  that  a  release  has  occurred,  but  without 
knowing  the  characteristics  of  the  source  (which  include  things  like  the  rate  and  duration  of 
the  release,  source  location,  etc.),  the  prediction  of  the  dispersion  of  the  contaminant  in  the 
atmosphere  cannot  be  made.  Indeed,  the  principal  impediment  and  critical  capability  gap 
in  current  emergency  management  efforts  directed  at  terrorist  incidents  involving  the  covert 
release  of  a  CBR  agent  into  the  atmosphere  is  the  determination  of  the  unknown  source 
characteristics  following  event  detection  by  a  network  of  CBR  sensors.  This  realization  has 
provided  the  impetus  for  our  recent  research  efforts  directed  towards  a  solution  of  the  source 
reconstruction  problem. 

Although  the  problem  of  the  forward  prediction  of  the  concentration  field  resulting  from 
the  dispersion  of  a  contaminant  released  in  a  complex  turbulent  flow  (with  a  known  source) 
has  received  considerable  attention  from  researchers,  the  “reverse”  prediction  of  the  source 
location  and  strength  using  a  finite  number  of  noisy  concentration  data  obtained  from  an 
array  of  sensors  has  received  considerably  less  attention,  although  the  importance  of  the 
solution  of  this  problem  for  a  number  of  practical  applications  is  obvious.  Early  work  on 
inverse  source  determination  focussed  on  the  problem  of  source  strength  estimation  whereby 
the  location  of  the  source  was  assumed  to  be  known  exactly  a  priori.  For  example,  Hanna  et 
al.  [1]  used  three  different  types  of  dispersion  models  to  estimate  the  emission  rate  (source 
strength)  for  the  localized  sources  used  in  the  Project  Prairie  Grass  experiments.  Gordon 
et  al.  [2]  and  Wilson  and  Shum  [3]  applied  a  forward-time  Lagrangian  stochastic  trajec¬ 
tory  model  in  conjunction  with  a  mass  balance  technique  to  estimate  the  rate  of  ammonia 
volatilization  from  field  plots  (viz.,  area  sources  on  the  ground  surface),  while  Flesch  et 
al.  [4]  used  a  backward-time  Lagrangian  stochastic  dispersion  model  to  estimate  the  emis¬ 
sion  rate  from  a  sustained  surface  area  source  over  simple  terrain.  Kaharabata  et  al.  [5] 
used  an  approximate  solution  to  the  advection-diffusion  equation  to  determine  the  source 
strength  from  microplots  over  an  open  field,  using  single  and  multipoint  measurements  of 
the  concentration  downwind  of  the  contaminated  field  plot.  For  longer-range  dispersion, 
Robertson  and  Langner  [6]  estimated  the  emission  rates  for  the  European  Tracer  Experi¬ 
ment  (Nodop  et  al.  [7])  using  a  variational  data  assimilation  procedure  to  minimize  a  cost 
function  defined  to  be  the  sum  of  squared  differences  between  the  model  solution  and  the 
concentration  measurements  over  a  24-h  assimilation  period  (window).  A  different  approach 
was  proposed  by  Siebert  and  Stohl  [8],  who  attempted  to  reconstruct  the  source  strength 
using  Tikhonov  regularization.  Finally,  Skiba  [9]  used  an  adjoint  transport  model  to  assess 
emission  rates  from  various  industrial  plants  in  Mexico,  whose  locations  were  assumed  to 
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be  known  a  priori. 


A  number  of  researchers  have  focussed  on  the  problem  of  recovering  the  location  of  an  un¬ 
known  source  on  the  basis  of  available  concentration  measurements,  given  that  the  source 
strength  is  known  a  priori.  Khapalov  [10]  formulated  the  source  localization  problem  as 
a  nonlinear  identification  problem  based  on  the  diffusion  system  modeling  involving  a  dis¬ 
tributed  system  of  parabolic  type.  Another  optimization  approach  was  followed  by  Alpay 
and  Shore  [11]  who  described  an  “intelligent”  brute- force  method  for  source  localization, 
which  is  based  on  the  simulation  of  every  possible  scenario,  with  the  likely  source  location 
being  that  which  minimizes  an  appropriate  norm  between  the  model  predictions  of  con¬ 
centration  for  that  source  location  and  the  measured  concentration  data.  Matthes  et  al. 
[12]  employed  a  two-step  approach  for  source  localization  for  the  simple  (albeit  unrealistic) 
case  of  an  emitted  substance  transported  by  homogeneous  advection  and  isotropic  diffusion. 
This  method  was  designed  to  overcome  convergence  difficulties  that  have  been  experienced 
by  various  researchers  using  gradient-based  optimization  methods  for  formulation  of  the 
source  localization  problem. 

A  simultaneous  estimation  of  the  source  location  and  strength  was  formulated  by  Alpay 
and  Shore  [11],  who  introduced  the  dual  system  method  whereby  the  forward  governing 
equation  for  the  concentration  is  replaced  with  the  adjoint  equation  for  an  adjoint  state 
variable.  A  similar  approach,  in  the  context  of  dispersion  on  a  global  scale,  was  followed 
by  Pudykiewicz  [13],  who  inferred  source  parameters  (location,  emission  rate,  and  time  of 
release)  by  using  the  influence  function  obtained  as  a  solution  of  the  adjoint  tracer  transport 
equation  to  estimate  the  source  location,  in  conjunction  with  the  Lagrange  duality  relation 
to  estimate  the  strength  of  the  source.  In  the  line  of  probabilistic  approaches,  a  Bayesian 
inferential  methodology  in  conjunction  with  the  adjoint  method  (in  which  the  transport  and 
diffusion  of  matter  in  turbulent  flows  was  solved  using  either  partial  or  stochastic  differential 
equations  backward  in  time  in  an  Eulerian  or  Lagrangian  description  of  turbulent  dispersion, 
respectively)  was  proposed  by  Yee  [14],  Yee  [15],  Keats  et  al.  [16],  and  Yee  [17]  for  the 
simultaneous  determination  of  all  source  parameters  for  dispersion  in  complex  environments 
(e.g.,  urban  dispersion,  long-range  dispersion,  etc.). 

In  all  previous  studies  cited  here,  with  the  exception  of  Yee  [17],  the  problem  of  identification 
of  source  parameters  was  restricted  to  one  source.  The  determination  of  the  characteristics 
of  multiple  localized  sources  was  briefly  addressed  in  Yee  [17],  but  the  approach  used  there 
assumed  that  the  number  of  sources  was  known  a  priori .  The  problem  of  reconstruction  of 
an  unknown  number  of  sources  using  a  finite  number  of  noisy  concentration  measurements 
obtained  from  an  array  of  sensors  is  of  great  interest  for  current  emergency  and  retrospective 
management  efforts  directed  at  terrorist  incidents  involving  the  release  of  CBR  agents  in 
a  densely  populated  urban  center.  In  this  application,  the  array  of  CBR  sensors  measures 
the  superposition  of  concentrations  arising  from  multiple  interfering  sources  in  which  both 
number  of  sources  and  the  characteristics  (e.g.,  location,  emission  rate)  of  each  source  are 
unknown  a  priori.  The  determination  of  the  number  of  sources  is  a  difficult  problem  and 
involves  some  form  of  model  selection.  In  this  paper,  we  formulate  a  Bayesian  inferential 
scheme  for  the  joint  determination  of  the  number  of  contaminant  sources  and  the  parame¬ 
ters  for  each  source,  given  a  finite  number  of  concentration  measurements  made  by  an  array 
of  sensors.  The  Bayesian  formulation  of  a  solution  for  this  problem  provides  a  computa¬ 
tional  challenge  owing  to  the  fact  that  the  resulting  integrals  of  the  posterior  distribution 
of  the  source  parameters  over  the  hypothesis  space  are  not  analytically  tractable,  and  stan- 
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dard  numerical  methods  for  integration  cannot  be  applied  to  give  accurate  results  due  to 
the  potential  high  dimensionality  of  the  hypothesis  space.  To  overcome  this  problem,  our 
proposed  procedure  uses  a  reversible-jump  Markov  chain  Monte  Carlo  method  to  estimate 
the  number  of  unknown  sources  and  the  relevant  parameters  (e.g.,  emission  rate,  source 
location,  etc.)  for  each  source. 


Concentration  data  model:  source-receptor  relationship 


Here  we  present  our  model  for  the  mean  concentration  “seen”  by  a  sensor  and  the  actual 
concentration  measurements  made  by  a  sensor,  the  latter  of  which  must  necessarily  include 
several  sources  of  “noise” . 


Model  for  mean  concentration 

To  solve  the  source  reconstruction  problem  using  Bayesian  probability  theory,  it  is  neces¬ 
sary  to  relate  the  hypotheses  of  interest  concerning  the  source  distribution  to  the  available 
concentration  data  measured  by  the  array  of  sensors.  This  is  a  source-receptor  relationship 
that  describes  the  transport  of  a  substance  (e.g.,  CBR  material)  through  the  atmosphere 
after  it  has  been  released  from  a  known  source  distribution  and,  as  such,  provides  a  pre¬ 
diction  of  the  average  concentration  of  the  substance  in  a  small  volume  centered  at  any 
given  spatial  location  during  any  given  time  interval.  With  reference  to  the  source-receptor 
relationship,  the  description  of  the  turbulent  transport  of  a  “marked”  fluid  element  (viz.,  a 
fluid  element  marked  by  the  contaminant  substance  as  it  moves  through  any  portion  of  the 
source  distribution)  in  a  Largrangian  reference  frame  is  perhaps  the  most  natural  and  sim¬ 
plest  description  of  this  physical  phenomenon.  This  involves  modeling  dispersion  through 
the  random  walks  for  “marked”  fluid  elements. 

Let  pi(x,  f|x',  t')(ix'  be  the  probability  that  a  fluid  element,  displaced  by  the  background 
flow  field  from  within  an  infinitessimal  volume  dx'  of  the  point  x'  (at  time  t'),  passes  through 
the  point  x  (at  time  t  where  t  >  t')  in  a  statistical  ensemble  of  turbulent  flows.  If  ^(x',  t') 
is  the  source  density  function  at  the  point  x'  (time  t'),  the  concentration  C  at  time  t  and 
location  x  of  the  contaminant  released  within  T)  (space  occupied  by  the  fluid)  up  to  this 
time  is  given  as  follows: 


C'(x,t)=  f  f  pi(x,  t|x',  t')5(x',  t')  dt' dx'.  (1) 

J  T>  J  —  oo 

The  mean  concentration  C{xr,tr)  “seen”  by  a  sensor  at  receptor  location  x^  and  time  R 
corresponds  to  an  average  of  C'(x,f)  over  the  spatial  volume  and  averaging  time  of  the 
sensor,  so 

C'(xr,tr)=  /  dt  dxC'(x,  t)/i(x,  t|Xr,  tr)  =  (C*, /l)(Xj.,  R).  (2) 

Jo  Jt> 

Here,  /i(x,  t|xr,R)  is  the  spatial-temporal  filtering  (or,  response)  function  (of  x  and  t)  for 
the  concentration  measurement  made  by  the  sensor  at  location  x^  and  at  time  R;  and, 
V  X  [0,  T]  is  the  space-time  volume  (of  the  sensor)  over  which  the  averaging  is  taken.  Note 
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that  in  Equation  (2),  we  expressed  the  mean  concentration  tr)  “seen”  by  a  sensor  as 

the  inner  (or  scalar)  product  {C,  h)  of  the  concentration  C  and  sensor  response  function  h. 

In  complex  turbulent  flows  (e.g.,  atmospheric  flow),  there  are  no  known  analytical  solu¬ 
tions  for  the  transition  probability  density  function  (PDF)  pi(x,  t|x',  t')  required  for  the 
determination  of  the  concentration  of  tracer  [cf.  Equation  (1)].  Because  the  exact  form  of 
Pi(x,  t|x',  t')  is  not  known,  we  need  to  construct  some  approximation  for  this  PDF.  In  the 
context  of  Lagrangian  stochastic  trajectory  simulation  methods,  the  transition  PDF  for  a 
fluid  element  at  position  x  at  time  t,  which  at  the  earlier  time  t'  was  at  position  x',  is  de¬ 
termined  as  the  solution  to  the  following  Ito  stochastic  differential  equation  (see,  Thomson 
[18]): 

dX(t)  =  U(t)  dt, 

^  /9  (3) 

dV{t)  =  a(X(t) ,  U(t) ,t)dt+  (Coe(X(t) ,  t)y  dW{t) , 

where  X(t)  =  {Xi{t))  =  X2{t),  X^it))  and  U(t)  =  {Ui{t))  =  {Ui{t),U2{t),U3{t)) 

are  the  (Lagrangian)  position  and  velocity,  respectively,  of  a  “marked”  fluid  element  (or, 
particle)  at  time  t  (marked  by  the  source  as  the  fluid  element  passes  through  it  at  some 
earlier  time  t'),  so  (X,U)  determines  the  state  of  the  fluid  particle  at  any  time  t  after 
its  intial  release  from  the  source  distribution  S.  In  Equation  (3),  Cq  is  the  Kolmogorov 
universal  constant  (associated  with  the  Kolmogorov  similarity  hypothesis  for  the  form  of  the 
Lagrangian  velocity  structure  function  in  the  inertial  subrange);  e  is  the  mean  dissipation 
rate  of  turbulence  kinetic  energy;  dW{t)  =  {dWi{t))  =  {dWi{t) ,  dW2{t) ,  dW3{t))  are  the 
increments  of  a  vector- valued  (three-dimensional)  Wiener  process  (i.e.,  they  have  a  Gaussian 
distribution  with  zero  mean  and  variance  dt  (infinitessimal  time  step) ,  and  non-overlapping 
increments  are  statistically  independent);  and  a  =  (a^)  =  (ai,  02,  03)  is  the  drift  coefficient 
vector. 


The  model  is  fixed  by  the  choice  of  the  drift  coefficients  (i  =  1,2,3).  The  most  important 
criterion  for  the  determination  of  is  the  well-mixed  criterion  [18]  which  states  that  an 
initially  well-mixed  distribution  of  particles  should  remain  so.  One  possible  model  for 
that  satisfies  the  well-mixed  criterion  was  given  by  Thomson  [18]  for  the  case  of  Gaussian 
turbulence  in  three  dimensions  in  which  the  (unconditional)  PDF  of  the  velocity  fluid 
elements  (both  marked  and  unmarked)  has  the  form: 


^e(u) 


1 

- T7? 

(27r)3/2(detrij) 


(4) 


where  Ui  is  the  mean  Eulerian  velocity  and  Tij  =  {ui  —  Ui){uj  —  U j)  is  the  Reynolds  stress 
tensor  for  the  underlying  turbulent  flow  (here,  tensor  notation  with  the  Einstein  summation 
convention  is  used  and  an  overline  is  used  to  denote  an  ensemble  average).  It  should  be 
noted  that  there  is  an  implicit  dependence  of  Pe  on  x  and  t  through  its  dependence  on 
Ui  =  Ui(x,t)  and  Tij  =  Tij(x,t)-  In  this  context,  Thomson’s  solution  for  ai  that  is 
consistent  with  the  well-mixed  criterion  assumes  the  following  form: 

ai  =  -^{Coe)T-f}{uk-Uk)  +  ^,  (5a) 
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where 


_  1  dVu  du^  -  du, 
Pe  2  dxi  dt  ^  dxi 


This  model  is  general  enough  to  allow  for  almost  any  flow  complexity. 

Equation  (3),  in  conjunction  with  Equations  (5a)  and  (5b),  is  a  forward-time  Lagrangian 
trajectory  simulation  model.  In  this  source-oriented  approach,  “marked”  fluid  elements 
with  initial  space-time  coordinates  (x',  t')  are  sampled  from  a  space-time  density  function 
proportional  to  the  (prescribed)  source  distribution  5(x',t').  The  forward  Lagrangian  tra¬ 
jectories  {t  >  t')  of  these  “marked”  fluid  elements,  which  emanate  from  the  source  and  move 
toward  the  receptor,  are  determined  in  accordance  to  Equation  (3).  This  algorithmic  proce¬ 
dure  permits  the  evaluation  of  the  mean  concentration  “seen”  by  a  sensor  at  any  arbitrary 
receptor  space-time  point  {'x^pr)  in  accordance  to  Equation  (2).  Note  that  evaluation  of 
the  inner  product  C{xr,tr)  =  {C,  h){xr,tr)  involves  two  basic  functions:  the  source  distri¬ 
bution  S'(x',t')  and  the  detector  space-time  filtering  function  h{x,t\xr,tr)-  Unfortunately, 
the  source-oriented  approach  described  here  is  too  computationally  expensive  for  use  in  a 
Bayesian  inferential  approach  for  multiple  source  reconstruction,  because  sampling  from  the 
posterior  distribution  for  the  source  parameters  (described  later  in  this  report)  requires  a 
large  number  of  forward  source-receptor  relationships  to  be  determined,  each  of  which  in¬ 
volves  the  solution  of  the  stochastic  differential  equation  given  by  Equation  (3).  Therefore, 
for  the  Bayesian  inversion  of  concentration  data  to  be  practical,  fast  and  efficient  techniques 
are  required  for  the  determination  of  the  source-receptor  relationship  within  the  context  of 
sampling  in  the  hypothesis  space  of  the  source  required  for  posterior  inference. 


In  view  of  this,  it  is  much  more  attractive  to  apply  a  receptor-oriented  approach  for  repre¬ 
sentation  of  the  source-receptor  relationship,  owing  to  the  fact  that  this  approach  provides 
a  fast  technique  for  predicting  mean  concentration  data  at  a  given  receptor  location  for 
arbitrary  hypotheses  about  the  source  distribution.  As  a  consequence,  we  adopt  and  de¬ 
velop  this  formulation  of  the  source-receptor  relationship,  as  it  provides  a  fast  and  reliable 
framework  for  Bayesian  inference  in  the  context  of  source  reconstruction.  To  this  purpose, 
we  consider  a  representation  for  the  source-receptor  relationship  that  is  dual  to  the  one 
given  by  Equation  (2)  as  follows: 


C(xr,tr)  =  /  dt  dxC(x,t^h{x,t\Xr,tr) 

Jo  Jt> 


tr-  (6) 

=  J  dt'  J  dx'C*  (x' ,t'\Xr,tr)S{x' ,t')  =  {C*  ,  S){Xr,tr), 

— oo  T) 


where  C* (x' ,t'\xr,tr)  is  an  adjunct  (dual)  concentration  at  space-time  point  {x',t')  asso¬ 
ciated  with  the  sensor  concentration  data  at  location  x^  and  time  p  (with  t'  <  p)-  We 
note  that  C*(x',  t'\xr,tr)  is  explictly  constructed  so  that  it  verifies  the  duality  relationship 
implied  by  Eq.  (6);  namely,  {C,  h)  =  (C*,  S). 
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A  backward-time  Lagrangian  trajectory  simulation  model,  that  is  dual  to  the  forward-time 
Lagrangian  trajectory  simulation  model  given  by  Equation  (1),  can  be  constructed  for  the 
computation  of  C*  so  that  it  exactly  satisfies  the  duality  relationship  {C,  h)  =  {C* ,  S).  To 
this  end,  suppose  the  backward-time  Lagrangian  trajectory  model  is  defined  as  the  solution 
to  the  following  stochastic  differential  equation: 

dX.\t')  =  lj\t')dt', 

(7) 

dV\t')  =  dt'  +  {Coe{X.\t'),t')y^^  dW{t'), 

where  t'  <  tr  and  at  any  given  time  t' ,  (X^,U^)  is  a  point  in  the  phase  space  along  the 
backward  trajectory  of  the  “marked”  fluid  element  (here  assumed  to  be  marked  or  tagged 
as  a  fluid  particle  which  at  time  passed  through  the  spatial  volume  of  the  detector 
at  location  x^).  The  displacement  statistics  of  “marked”  fluid  elements  released  from  the 
receptor  location  at  time  tr  can  be  used  to  compute  C'*(x',  t'\xr,tr)  (which  is  interpreted 
here  as  a  function  of  x'  and  t').  It  can  be  shown  (Thomson  [18],  Flesch  et  ah  [4])  that  C* 
obtained  from  Equation  (7)  for  a  detector  with  the  filtering  function  h  and  C  obtained  from 
Equation  (2)  for  a  release  from  the  source  density  S  is  exactly  consistent  with  the  duality 
relationship  (C,  h)  =  (C*,  S)  if  and  only  if  in  Equation  (7)  is  related  to  a  in  Equation  (3) 
through  the  following  gauge  transformation: 


d 

a!’{x,u,t)  =  ai(x,u,t)  -  C'oe(x,  t)— InPB(u) 

OUi 


(8) 


This  form  of  the  drift  coefficient  for  the  backward-time  Lagrangian  stochastic  model  is  a 
consequence  of  Thomson’s  well-mixed  criterion  and  guarantees  that  the  duality  relationship 
between  C*  and  C,  for  any  source  distribution  S  and  detector  response  function  /i,  is  exactly 
satisfied.  Substituting  Equations  (4)  and  (5a)  in  Equation  (8),  an  explicit  form  for  a(’  can 
be  obtained  as 

a,^  =  i(Coe)r-;(u,-I7,)  +  |l.  (9) 

Note  that  a(’  differs  from  by  virtue  of  a  sign  change  on  the  first  term  (or,  “fading”  memory 
term)  on  the  right-hand-side  of  Equation  (9).  The  second  term  (“drift”  term  due  to  spatial 
inhomogeneity  and/or  temporal  non-stationarity  in  the  velocity  statistics  of  the  turbulent 
flow)  on  the  right  remains  invariant  under  a  time  reversal  operation.  This  suggests  that  the 
“fading”  memory  and  “drift”  terms  transform  differently  under  a  time  reversal  operation. 


Model  for  concentration  observations 


The  models  described  above  provide  predictions  for  the  “ideal”  mean  concentration  seen  by  a 
sensor  at  the  receptor  space-time  point  (x^,  tr)-  The  actual  concentration  data  measured  by 
the  sensor  will  not  usually  agree  with  the  concentration  predicted  by  the  model  owing  to  the 
noise  process  imposed  on  the  concentration  data,  which  by  its  very  nature  is  expected  to  have 
a  very  complicated  structure.  To  this  purpose,  it  is  assumed  that  the  actual  concentration 
data  available  from  the  sensor  array  were  measured  at  a  finite  number  of  sensor  locations  and 
at  a  finite  number  of  time  points  at  each  sensor  location.  The  actual  concentration  datum 
dj  j(»)  acquired  by  the  sensor  at  receptor  location  x^.  and  at  time  (i  =  1, 2, . . . ,  and 
j  =  1,2,...,  ,  where  is  the  number  of  sensors  and  '  is  the  number  of  time  samples 
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measured  at  the  i-th  sensor)  is  assumed  to  be  the  sum  of  a  modeled  mean  concentration 
signal  C(xj. . ,  ;  0)  and  “noise”  so 


(0 


=  ci^r^4\e)  + 


(0  ) 


(10) 


where  0  is  an  appropriate  parameter  vector  describing  the  source  distribution  5;  and, 
C{xr,tr',Q)  is  the  modeled  mean  concentration  at  location  x,.  and  time  tr,  determined  in 
accordance  to  Equation  (6)  for  a  source  distribution  characterized  by  parameter  vector  0. 
For  simplicity  of  notation,  the  variables  in  Equation  (10)  which  are  indexed  or  labelled  by 
(i,  will  be  ordered  in  some  regular  and  convenient  manner  (e.g.,  lexicographic  ordering) 
and  this  collection  will  be  indexed  by  J  ( J  =  1,  2, . . . ,  A^,  with  N  =  being  the 

total  number  of  measured  concentration  data).  Then,  we  can  write  the  observational  model 
as  follows: 

dj  =  Cj{Q)  +  ej,  J=1,2,...,N,  (11) 

where  Cj(0)  =  C(x^^,  0). 


In  Equation  (11),  the  ej  is  a  noise  term  representing  the  uncertainty  in  dj.  In  general,  ej 
consist  of  errors  (e.g,  input,  stochastic,  and  measurement)  and  any  real  signal  in  the  data 
that  cannot  be  explained  by  the  model.  The  random  error  ej  can  be  split  into  four  terms 
as  discussed  by  Rao  [19],  so 


ej  =  +  PP  +  pp  +  pp. 


(12) 


The  first  term  pp  of  the  error  corresponds  to  model  error  arising  from  uncertainties  in  the 
representation  of  various  physical  processes  in  the  dispersion  model  used  to  predict  the  mean 
concentration.  The  second  term  rjj  ^  describes  the  input  error  arising  from  uncertainties 
in  the  values  of  empirical  parameters  and/or  specification  of  the  input  meteorology  (initial 
and  boundary  conditions)  used  by  the  dispersion  model.  The  third  term  pp  of  the  error  is 
the  stochastic  uncertainty  arising  from  the  turbulent  nature  of  the  atmosphere,  which  gives 
rise  naturally  to  random  concentration  fluctuations  in  hazardous  gas  releases.  Finally,  pp 
describes  the  noise  inherent  in  the  sensor  (essentially  measurement  or  instrument  error) . 


Rao  [19]  discusses  the  nature  of  these  four  types  of  error  with  respect  to  characterization 
of  uncertainties  in  atmospheric  dispersion  models,  and  provides  a  comprehensive  review  of 
sensitivity  and/or  uncertainty  analysis  methods  that  have  been  used  to  quantify  and  reduce 
them.  In  this  paper,  all  the  various  error  contributions  to  the  noise  term  are  simply  lumped 
together  and  denoted  by  ej  [see  Equations  (11)  and  (12)].  It  is  assumed  that  the  observer 
does  not  have  a  detailed  knowledge  of  the  probability  distribution  of  the  noise  (aggregate 
error),  other  than  that  the  observer  has  an  estimate  for  the  expected  scale  of  variation  of 
the  noise.  More  specifically,  it  is  assumed  that  the  noise  scale  parameter  associated  with  ej 
(aggregate  error  of  the  J-th  concentration  observation)  is  provided  in  the  form  of  a  (finite) 
variance  Uj. 

In  this  paper,  it  is  considered  that  the  model  concentration  Cj{Q)  in  Equation  (11)  results 
from  the  release  of  a  contaminant  (e.g.,  CBR  material)  from  Ng  transient  point  sources.  It  is 
assumed  that  the  k-th  source  was  located  at  x^  ^  and  that  this  source  was  activated  (turned 
on)  at  time  and  deactivated  (turned  off)  at  a  later  time  Tp  Furthermore,  over  the  time 
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period  from  to  T^,  the  source  was  assumed  to  be  releasing  contaminant  at  a  constant 
emission  rate  Qk  {k  =  1,2, ...  Ng).  Consequently,  the  (unknown)  source  distribution  is 
assumed  a  priori  to  have  the  following  form: 


5(x,  t)  =  Y,  QkS{^  -  ^■s,k)  U{t  -  T^)  -  U{t  -  T^) 


k=l 


(13) 


where  5{-)  and  [/(•)  are  the  Dirac  delta  and  Heaviside  unit  step  functions,  respectively.  It  is 
convenient  to  define  a  source  parameter  vector  6 corresponding  to  the  source  distribution 
S  of  Equation  (13)  as  On,  =  {xg^i,Tl  ,Tg  ,  Qi, . . .  ,Xg^N,,T^%T^%  QnJ  G  Further¬ 

more,  let  0  =  {Ng,6N,).  If  we  substitute  Equation  (13)  into  Equation  (6),  the  model 
concentration  C  j{Q)  “seen”  by  the  sensor  at  location  x^j  and  time  trj  is  given  explicitly 

by 


Cj{Q)  =  'S^Qk  /  C*{yis^k,ts\y.rj,trj)dts. 

rrt  JTt 


(14) 


The  problem  of  source  reconstruction  reduces  to  the  following:  given  the  observed  vector 
of  concentration  data  D  =  {di,  0,2, . . . ,  cIn),  the  objective  is  to  estimate  Ng  and  On,  or, 
equivalently,  0. 


Bayesian  analysis 


To  determine  the  number  of  sources  and  reconstruct  the  characteristics  (e.g.,  location,  emis¬ 
sion  rate)  of  each  source  (parameters  encapsulated  in  0),  given  the  vector  of  concentration 
observations  D  obtained  from  the  array  of  sensors,  the  techniques  of  Bayesian  inference  are 
employed.  The  Bayesian  framework  is  attractive  in  the  current  context  because  it  provides 
a  rigorous  mathematical  foundation  for  making  inferences  about  the  source  parameters 
and,  as  a  consequence,  provides  a  rigorous  basis  for  quantifying  the  uncertainties  in  the 
estimated  source  parameters.  Bayesian  inference  can  be  obtained  from  the  product  rule  of 
probability  calculus,  the  latter  of  which  can  be  derived  rigorously  starting  with  the  formu¬ 
lation  of  a  small  number  of  desiderata  required  to  define  a  rational  theory  of  inference  as 
first  enunciated  by  Cox  [20],  with  a  more  complete  treatment  given  by  Jaynes  [21]  in  his 
definitive  treatise.  This  formulation  leads  to  the  ordinary  rules  (sum  and  product  rules)  of 
probability  calculus  and  implies  that  every  allowed  (consistent)  theory  for  inference  must 
be  mathematically  equivalent  to  probability  theory,  or  else  inconsistent  (no  other  calculus 
is  admissible  for  inference  that  is  consistent  with  the  above  mentioned  desiderata). 

The  basic  relationship  quantifying  parameter  inference  is  Bayes’  rule,  which  in  the  context 
of  the  current  problem  can  be  expressed  as  follows: 

p(0iD,  i)p{B\i) = p(D|0,  mm-  (15) 

The  factors  in  Equation  (15)  are:  (1)  p{Q\I)  is  the  prior  probability  density  function  (PDF) 
of  the  source  parameter  vector  0  that  encapsulates  our  state  of  knowledge  of  the  param¬ 
eters  before  the  receipt  of  the  concentration  measurements;  (2)  p(Dj0,/)  is  termed  the 
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likelihood  function  when  considered  as  a  function  of  0,  but  is  known  as  the  sampling  distri¬ 
bution  when  considered  as  a  function  of  D  (the  latter  being  the  probability  of  observing  the 
concentration  data  D  when  given  the  source  parameters  0);  (3)  p(D|/)  is  termed  the  evi¬ 
dence  (also,  frequently  referred  to  as  the  prior  predictive  or  the  marginal  likelihood)  ;  and, 
(4)  p(0|D,/)  is  the  posterior  probability  density  function  of  the  parameters  0  of  interest, 
that  corresponds  to  the  update  of  p{Q\I)  incorporating  the  knowledge  gained  about  0  after 
receipt  of  the  concentration  observations  D.  All  the  terms  here  are  to  be  interpreted  given 
the  background  (contextual)  information  I  (e.g.,  background  meteorology,  source-receptor 
relationship,  etc.)  that  dehnes  the  source  reconstruction  problem.  This  is  the  meaning  of  I 
after  the  vertical  bar  “|”  which  is  used  here  to  denote  “conditional  upon”. 

Using  Bayesian  inference,  all  information  in  the  measured  concentration  data  relevant  to 
the  problem  of  estimating  the  source  parameter  vector  0  is  summarized  in  the  posterior 
PDF  p(0|D,  I),  which  by  virtue  of  Equation  (15),  is  given  by 

p(0|D,/)ocp(D|0,/)p(0|/).  (16) 

It  is  noted  that  for  parameter  estimation,  the  evidence  p(D|/)  is  0-independent  and  simply 
plays  the  role  of  a  normalization  factor;  hence,  p(0|D,/)  is  expressed  simply  as  a  propor¬ 
tionality  in  Equation  (16).  To  proceed  further  in  the  specihcation  of  the  posterior  PDE,  we 
now  need  to  assign  functional  forms  for  the  prior  PDE  p(Q\I)  and  the  likelihood  function 
p(D|0,/). 

Eirst,  let  us  consider  the  assignment  of  the  prior  PDE  p{Q\I)  =  p{{Ng,9]\[J\I).  To  this 
purpose,  we  assume  the  logical  independence  of  the  various  source  parameters;  that  is  to 
say,  knowing  the  value  of  the  number  of  sources  tells  us  nothing  about  the  various  source 
characteristics  (e.g.,  location,  emission  rate),  knowing  the  characteristics  of  one  source  tells 
us  nothing  about  those  of  another  source,  knowing  the  location  of  a  particular  source  tells  us 
nothing  about  the  emission  rate  or  when  the  source  was  turned  on/off,  etc.  In  consequence, 
the  prior  probability  p{{Ns,  6]y^)\I)  factorizes  (by  repeated  application  of  the  product  rule 
of  probability  calculus)  as  follows: 

Ns 

p{{Ng,eNs)\I)=p{Ng\I)  ]Jp{Qk\I)p{^s,k\I)p{Tt\I)p{T^\Ttl).  (17) 

k=l 


We  now  discuss  the  assignment  of  each  of  the  (component)  prior  distributions  in  Equa¬ 
tion  (17).  The  prior  distribution  on  Ng  (number  of  sources)  is  chosen  to  be  a  truncated 
Poisson  distribution  with  parameter  P,  where  P  is  the  expected  number  of  sources,  so 


p{Ng\I)  oc  P'^^ 


exp(-P) 

Ng\ 


(18) 


where  A^s,max  is  the  maximum  number  of  sources  and  ]Ia{x)  is  the  indicator  function  for 
set  A  (viz.,  Ia{x)  =  1  if  x  G  A,  and  Ma{x)  =  0  if  x  ^  A).  Here  and  elsewhere,  irrelevant 
normalization  constants  in  functional  forms  of  probability  distributions  will  be  ignored. 
Given  Ng,  the  prior  distributions  for  Qk,  and  {k  =  1,2,...,W)  are  each 

chosen  to  be  uniform  over  an  appropriate  domain  of  dehnition.  More  specihcally,  the  prior 
distribution  for  is  uniform: 


p{Qk\I)  =  ^(0,Qn,s..  ){Qk)/ Q  max;  k  —  1,2,...,  Ng , 


(19) 
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where  Qmax  is  an  upper  bound  for  the  expected  emission  rate.  The  prior  distribution  for 
is  chosen  to  be  uniform  over  some  large  region  D  C 

=  ^vi^s,k)/V(T>),  k  =  1,2, . .  .,Ns,  (20) 

where  ViV)  is  the  volume  of  the  region  T>.  Finally,  the  prior  distributions  for  and  Tg 
are  assigned  uniform  distributions  with  the  following  forms: 

p{T^\I)  =  l(0,T_)(r^)/r,,ax,  k  =  l,2,...,N,-,  (21) 

and 

p(rg"|r^, I)  oc  -t^),  k  =  i,2,...,N,.  (22) 

Here,  is  the  upper  bound  on  the  time  at  which  the  source  was  turned  on  and  turned 
off.  Note  that  the  time  that  the  source  is  turned  off  must  necessarily  occur  after  it  has 
been  turned  on,  and  this  information  is  encoded  in  form  of  the  prior  distribution  for 
given  by  Equation  (22),  where  the  distribution  is  seen  to  be  conditioned  on  .  In  view  of 
Equations  (18)  to  (22),  the  prior  distribution  for  0  given  by  Equation  (17)  now  assumes  the 
following  explicit  form  (up  to  a  normalizing  constant  that  does  not  depend  on  the  source 
parameters): 


Ns 

k=l 

(23) 

(2^max  —  ) 

Next,  consider  the  assignment  of  a  functional  form  for  the  likelihood  function  p(D|0,/). 
The  likelihood  function  is  equivalent  to  the  direct  probability  for  the  concentration  data 
D  given  the  source  parameters  0.  In  the  absence  of  a  detailed  knowledge  of  the  noise 
distribution  ej  [cf.  Equation  (11)],  other  than  that  it  has  a  finite  variance  iTj,  the  application 
of  the  principle  of  maximum  entropy  (Jaynes  [21])  informs  us  that  a  Gaussian  distribution 
is  the  most  conservative  choice  for  the  direct  probability  of  the  data  D  (or,  equivalently, 
for  the  noise  e  =  (ci,  62, . . . ,  cat)).  It  is  important  to  note  that  the  entropy  of  the  noise 
PDE  (which  is  equivalent  to  the  likelihood  function)  is  a  measure  of  the  size  of  the  basic 
support  set  of  the  distribution  or  ‘volume’  occupied  by  the  sensibly  probable  noise  values. 
Choosing  a  distribution  for  the  noise  that  provides  the  largest  support  set  permitted  by 
the  available  information  allows  the  largest  range  of  possible  variations  in  the  noise  values 
consistent  with  that  information,  implying  the  most  conservative  estimates  for  these  values. 
Assigning  a  Gaussian  distribution  for  the  noise  using  the  maximum  entropy  principle  makes 
no  statement  about  the  true  (unknown)  sampling  distribution  of  the  noise  (which  has  been 
indicated  earlier  to  have  a  highly  complicated  structure  for  the  current  problem).  Rather,  it 
simply  represents  a  maximally  uninformative  state  of  knowledge,  a  state  of  knowledge  that 
reflects  what  the  observer  knows  about  the  true  noise  in  the  data  (namely,  the  mean  and 
variance  of  the  noise,  with  all  other  properties  of  the  noise  being  irrelevant  to  the  inference 
since  these  are  unknown  to  the  observer).  Erom  these  considerations,  the  likelihood  function 
for  our  problem  has  the  following  form  [in  light  of  Equation  (11)]: 


10 
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where 


N 


X^{&)  =  E 


J=1 


dj-Cj{e\j 


o-j 


(246) 


Inserting  Equations  (23)  and  (24)  into  Equation  (16),  the  posterior  distribution  for  0  has 
the  following  form,  where  all  superfluous  constants  that  are  independent  of  Ng  or  have 
been  absorbed  into  the  constant  of  proportionality  (or,  normalization  constant): 


p(0|D,/)oc 


exp 


j=i 


Uj=i 

. 


0-j 


Ngl 


N, 


k=l  f-^max  J^b) 


recalling  that  the  source  parameter  vector  0  =  {Ng,  ^nJ-  Here,  C'j(0)  is  determined  in  ac¬ 
cordance  to  Equation  (14)  with  the  adjunct  “concentration”  field  C'*(x',  t'|x^,  t^)  predicted 
using  a  backward-time  Lagrangian  stochastic  model  [cf.  Equations  (7),  (9),  and  (5b)].  The 
posterior  distribution  for  0  provides  the  full  solution  for  the  multiple  source  reconstruction 
problem.  Inferences  on  the  values  of  the  source  parameters  are  based  on  upon  this  posterior 
distribution.  This  posterior  distribution  may  be  summarized  by  various  statistics  of  interest 
such  as  the  posterior  mean  of  each  source  parameter,  say  9]^  (which  can  be  an  emission 
rate  Qk,  or  source  location  x^^fc,  or  time  at  which  the  source  was  turned  on  (off)  (Tg  ) 
for  the  k-th.  source  for  a  source  distribution  consisting  of  Ng  sources  with  k  =  \,  2, Ng): 


=£:[6ivjD]  =  J  e^^^p{{Ng,eNj\T>,I)deN^,  (26a) 

where  £[■]  denotes  mathematical  expectation.  A  measure  of  uncertainty  of  this  estimate  of 
is  the  posterior  standard  deviation 

^(^k)  =  =  (I (0k  '  .  (266) 

Alternatively,  a  p%  credible  (or,  highest  posterior  distribution)  interval  that  contains  the 
source  parameter  0*  with  p%  probability,  with  the  lower  and  upper  bounds  of  the  interval 
specified  such  that  the  probability  density  within  the  interval  is  everywhere  larger  than  that 
outside  it,  can  be  used  as  a  measure  of  the  uncertainty  in  the  determination  of  0*.  Einally, 
an  estimate  for  the  number  of  sources  can  be  obtained  from  the  maximum  a  posteriori 
estimate  as  follows: 

Ng  =  argmax  p{Ng\D,I),  (27) 

^  s 

where  p{Ng\'D,  I)  is  the  posterior  probability  of  the  number  of  sources,  given  the  concen¬ 
tration  data  and  the  contextual  background  information.  The  posterior  probability  for  the 
number  of  sources  is  a  marginal  posterior  probability,  where  all  the  parameters  we  are  not 
interested  in  (e.g.,  0Ar^)  are  removed,  through  the  process  of  integrating  over  these  param¬ 
eters  (in  a  process  called  marginalization). 
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Bayesian  computation  and  Markov  chains 


A  perusal  of  Equation  (25)  shows  that  the  posterior  distribution  is  highly  nonlinear  in  the 
source  parameters  associated  with  location  (x^  fc)  and  with  the  times  the  source  was  turned 
on  {T^)  or  off  (T^  ),  and  that  explicit  evaluations  of  Equations  (26)  and  (27)  (and  similar 
integrals  or  functionals  arising  in  Bayesian  analysis)  are  impossible.  In  view  of  this,  we  apply 
posterior  sampling  for  evaluation  of  these  integrals,  which  is  implemented  using  a  Markov 
chain  Monte  Carlo  (MCMC)  algorithm  (Wilks  et  al.  [22],  Gelman  et  al.  [23]).  A  MCMC 
algorithm  can  be  used  to  generate  samples  of  source  distribution  models  (characterized  by 
0)  from  the  posterior  distribution  in  Equation  (25).  Towards  this  objective,  a  Markov  chain 
)}  is  constructed  whose  stationary  (or,  invariant)  distribution  is  the 
posterior  distribution  p(0jD,  I)  of  the  parameters  0  =  {Ng,  On^)- 


All  quantities  of  interest,  such  as  posterior  means  of  various  source  parameters  and  various 
marginal  posterior  distributions,  can  be  estimated  by  sample  path  averages  of  the  Markov 
process  {0*^*^}.  Eor  example,  the  marginal  posterior  distribution  p(A^sjD,/),  which  is  re¬ 
quired  in  Equation  (27)  for  the  inference  of  Ng,  can  be  estimated  as  follows: 


p(iV,lD,/)  =  ^#{t:ArW=Ar4 


1 


t=i 


(28) 


where  is  the  number  of  samples  0*^^\  0*-^^  . . .  drawn  from  a  sampled  realization  of  the 
Markov  chain  and  J|jv,}(')  is  used  to  select  those  samples  {0i*i,  t  =  1,2, . . . ,  N^}  drawn  from 
the  Markov  chain  that  correspond  to  source  distribution  models  with  exactly  Ng  sources 
(viz.,  those  samples  with  Ng*^  =  Ng).  Similarly,  from  the  samples  generated  from  such  a 
Markov  chain,  the  posterior  means  of  various  source  parameters  9)^  [see  Equation  (26a)] 
can  be  estimated  by 


(29) 


The  difficulty  in  the  construction  of  a  Markov  chain  for  sampling  from  the  posterior  dis¬ 
tribution  in  Equation  (25)  resides  in  the  fact  that  Ng  is  an  unknown,  requiring  the  design 
of  a  Markov  chain  that  can  simultaneously  explore  both  parameter  and  model  spaces.  In 
the  current  problem,  we  need  to  consider  a  set  of  models  Ai  =  indexed  by 

the  integer  parameter  Ng  £  Mg  =  {1,2,...,  A^s,max})  each  with  a  source  parameter  vector 
Ons  ^  Here,  model  corresponds  to  all  permissible  source  distributions  having 

exactly  Ng  sources.  The  simultaneous  exploration  of  a  number  of  candidate  source  models 
involving  differing  numbers  of  sources  can  be  realized  using  the  reversible-jump  MCMC 
algorithm  that  was  originally  introduced  by  Green  [24]  in  a  Bayesian  model  determination 
setting. 

The  remainder  of  this  section  will  focus  on  the  design  of  a  reversible-jump  MCMC  algorithm 
for  multiple  source  reconstruction,  in  the  case  where  the  number  of  sources  is  unknown  a 
priori.  The  MCMC  sampling  required  to  address  this  problem  involves  three  basic  compo¬ 
nents:  (1)  within-model  (or,  propagation)  moves  in  the  Markov  chain  involving  updating 
the  state  of  the  chain  on  a  hypothesis  (source  parameter)  space  of  fixed  dimension  (Ng  fixed) 
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using  a  standard  Metropolis-Hastings  update;  (2)  application  of  a  reversible-jump  MCMC 
for  between-model  moves  in  the  Markov  chain  involving  “jumps”  (creation  or  annihilation  of 
sources)  in  the  state  of  the  chain  between  hypothesis  (source  parameter)  spaces  of  different 
dimensions  corresponding  to  source  distributions  having  different  numbers  of  sources  {Ng 
variable);  and,  (3)  incorporation  of  parallel  tempering  and  Metropolis-coupled  MCMC  to 
increase  the  “mixing  rate”  of  states  of  the  Markov  chain  and  to  achieve  improved  compu¬ 
tational  efficiency.  We  will  now  describe  each  of  these  components  of  the  Markov  chain  in 
the  remainder  of  this  section. 


Propagation  moves  and  MCMC 


The  objective  of  MCMC  is  to  construct  a  Markov  chain  whose  stationary  distribution  is 
one  that  coincides  exactly  with  the  target  probability  density  function  that  we  are  trying  to 
sample  from  [e.g.,  in  our  case  that  target  probability  distribution  is  the  posterior  distribution 
p(0|D,  I)  given  in  Equation  (25)].  In  this  subsection,  we  consider  the  problem  of  design  of 
a  MCMC  algorithm  that  samples  from  p(0|D,  I)  for  fixed  Ng.  The  basic  MCMC  algorithm 
consists  of  two  components:  (1)  a  transition  or  proposal  distribution  function  (or  transition 
kernel)  T(0'|0);  and,  (2)  an  acceptance  probability  q;(0,  0').  These  two  components  are 
related  as  follows:  given  a  chain  in  the  current  state  0^*1  =  0  at  iteration  t,  a  proposed 
new  state  0'  =  ^(0)  =  V{{Ng,  On,))  =  {Ng,  is  drawn  from  some  proposal  (transition) 
distribution  T(0'|0)  and  this  new  point  is  accepted  as  the  new  state  of  the  chain  at  iteration 
{t  +  1)  with  the  standard  Metropolis-Hastings  (M-H)  acceptance  probability  (Gelman  et  al. 
[23])  given  by 


a(0,  0')  =  min 


’p(0lD,/)r(010O  j  ■ 


(30) 


Note  that  'P(-)  is  a  pure  propagation  operation  (move)  that  “translates”  9^^  G  to 

G  ,  while  keeping  Ng  fixed  (update  move  in  a  fixed  dimensional  hypothesis  space) . 
If  the  proposal  for  this  propagation  move  is  accepted,  then  the  new  state  at  iteration  (t-|- 1) 
is  00+^)  =  0';  otherwise,  00+^1  =  0.  Finally,  it  is  noted  that  the  M-H  algorithm  does  not 
require  the  normalization  of  p(0jD,  I)  be  known  [cf.  Equation  (30)]. 


To  proceed  further,  it  is  necessary  to  specify  appropriate  forms  for  the  proposal  density 
T(0'j0).  For  the  propagation  move,  Ng  is  fixed  and  the  update  occurs  only  for  9^^  ■  To  this 
purpose,  it  is  useful  to  partition  9]\f^  into  two  blocks  of  parameters  as  follows:  9n^  =  {9^,  9“^) 
where  9^  =  {Qi,Q2,  •  •  • ,  QivJ  G  and  9^  =  (x,,i,ri,ri, . . . ,  G 

This  particular  partitioning  of  the  source  parameter  vector  distinguishes  parameters  that 
are  related  linearly  (9^)  and  nonlinearly  (0^)  to  the  model  concentration  data  as  specified 
in  Equation  (14).  With  this  partitioning  of  the  source  parameter  vector,  a  cycle  of  two 
different  types  of  “blocked”  component  updates  (Gibbs  step  and  M-H  step)  involving  9^ 
and  9^  (respectively)  is  combined  in  sequence  to  form  a  single  iteration  of  the  MCMC 
sampler  for  generation  of  the  updated  state  of  the  Markov  chain. 


First,  let  us  consider  updating  the  emission  rates  Qk  of  the  various  sources  {k  =  1,2, Ng) 
which  we  have  collected  together  in  9^.  For  the  first  part  of  the  iteration  {t  +  1),  we  fix 
the  source  parameters  in  9^  =  (e.g.,  source  locations,  activation  times,  deactivation 

times)  to  the  values  obtained  in  the  previous  iteration  t  and  focus  on  the  re-sampling  of 
Qk  {k  =  1,2,...,W)-  We  choose  here  to  sample  the  emission  rates  one-at-time  using 
Gibbs  sampling.  The  Gibbs  sampler  updates  Qk  as  a  direct  draw  from  the  univariate 
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full  conditional  posterior  distribution  ,  D, /),  where  is  used  to  denote 

a  subvector  containing  all  the  components  of  9^  with  the  exception  of  Qk-  The  Gibbs 
sampler  is  simply  a  special  case  of  the  M-H  sampler  with  the  proposal  density  T(0'|0)  = 
p{Qk\9^i^,0‘^^*\H,  I)  [cf.  Equation  (30)]  for  the  update  of  the  Qk  component  of  0.  With 
this  proposal  density,  it  is  straightforward  to  show  that  the  M-H  acceptance  probability  for 
the  Gibbs  transition  is  exactly  unity  (all  draws  from  the  proposal  density  are  accepted). 
More  specifically,  the  re-sampling  of  Qk  {k  =  1,2, . . .,  W)  for  iteration  {t  +  1)  proceeds 
using  a  systematic  sweep  Gibbs  sampling  strategy  as  follows: 

simulate  ^  p{Qi\Q2  \  ■  ■  •  >  T), /); 

simulate  ~  p{Q2\Qf^"\  •  •  • ,  D,  /); 


simulateQS^+')  02  W,D,/);  (31) 

which  yields  0^1*+^)  =  {Qi~^^\Q2~^^\...,Q^^'^^'')  after  Ng  cycles.  It  should  be  noted 
that  Gibbs  sampling  can  be  applied  to  re-sample  Qk  because  the  univariate  conditional 
distribution  p{Qk\9^j,,  ,  I)  has  a  particularly  simple  form  from  which  it  is  possible 

to  obtain  Qk  as  a  direct  draw.  Indeed,  a  perusal  of  Equations  (25)  and  (14)  shows  that 
p{Qk\9^k^  T),  /)  is  simply  a  truncated  Gaussian  distribution  with  the  following  form: 


p(gfc|6']_fc,02W,D,/)  oc  exp  -Qk[J2  '^Qrfr  i^'^)  -  dj 


N 


Ns 


J=1  '-r=l 
r^k 

N 


(32a) 


j=i 


where  for  simplicity  of  notation  we  have  defined 


^mi 


C  (x^ /j,  jxj. j,  g j)  (its 


(326) 


Eollowing  the  update  of  Qk  for  iteration  (t  +  1)  in  accordance  to  Equation  (31),  the  second 
part  of  this  iteration  involves  updating  the  nonlinear  source  parameters  collected  in  0^  from 
021*1  to  0^1*+^).  To  accomplish  this  part  of  the  iteration,  the  emission  rate  parameters 
are  assumed  to  be  fixed  at  0^  =  0^1*+^!.  Eor  the  updating  of  these  parameters,  it  is  no 
longer  possible  to  use  a  Gibbs  sampler  because  the  resulting  univariate  full  conditional 
distribution  of  these  parameters  cannot  be  sampled  from  directly.  In  view  of  this,  we 
need  to  perform  a  M-H  step  to  update  these  parameters  with  an  appropriately  selected 
proposal  distribution.  Towards  this  objective,  we  sample  each  of  the  nonlinear  parameters 
(locations,  activation  times,  and  deactivation  times  of  the  various  sources)  using  a  M-H  step 
with  proposal  distribution  T{9‘^  |0^)  given  as  follows: 


et\elr^N{el(3l)  = 


[€  -  01] 


V^^k 


exp 


212 


k  =  l,2, 


ANg 


(33) 
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for  valid  0^  (e.g.,  if  the  prior  information  encoded  in  the  prior  distribution  requires  the 
parameter  to  lie  in  a  certain  interval,  then  a  proposal  for  this  parameter  that  lies  outside 
this  interval  is  rejected).  In  Equation  (33),  0^  {k  =  1,2, . .  .,4:Ns)  denotes  the  k-th.  compo¬ 
nent  of  9^.  The  proposal  distribution  here  provides  a  candidate  0^  that  is  a  perturbation 
of  the  current  value  of  0^,  which  is  obtained  by  drawing  from  a  Gaussian  distribution  with 
variance  .  Each  f3k  determines  the  scale  of  the  proposal  steps  for  the  nonlinear  parameter 
9l,  and  this  scale  needs  to  be  carefully  selected  in  order  to  ensure  an  appropriate  acceptance 
probability  for  the  proposed  (propagation)  move.  In  particular,  too  small  a  value  for  (3^  will 
lead  to  small  propagation  steps  that  are  almost  all  accepted  (but  result  in  a  very  slow  explo¬ 
ration  of  the  hypothesis  space) ,  whereas  too  large  a  value  for  f3k  will  result  in  an  excessively 
high  rejection  rate  for  the  proposed  move  (and,  in  so  doing,  also  severely  restrict  movement 
in  the  hypothesis  space) .  It  is  a  time-consuming  process  to  find  an  appropriate  scale  for  the 
proposal  distribution  for  each  parameter.  In  an  attempt  to  circumvent  this  problem,  we  use 
proposal  distributions  that  are  mixtures  of  four  Gaussian  distributions  centered  on  O^,  each 
with  the  form  given  by  Equation  (33).  However,  each  Gaussian  distribution  of  this  mixture 
has  a  different  standard  deviation  (3,  with  the  standard  deviations  chosen  in  such  a  manner 
as  to  cover  several  different  orders  of  magnitude.  This  strategy  ensures  that  a  near  optimal 
scale  for  the  step  size  is  used  at  least  some  of  the  time.  Because  the  proposal  distribution 
in  Equation  (33)  is  symmetric  (for  transition  between  Of,  and  Of  ),  the  acceptance  probabil¬ 
ity  in  Equation  (30)  for  the  transition  Of  —>■  Of  ,  with  the  posterior  probability  defined  in 
Equation  (25),  is  given  simply  by 

a(0,  0')  =  min  1 1,  exp  (^x^(0')  -  X^(0))^  |  ,  (34) 

where  0'  here  involves  only  the  change  of  the  k-th  component  in  the  subvector  0^. 

Equations  (31)  and  (33)  completely  specify  one  iteration  involving  the  sampling  of  a  new 
candidate  state  0]\f^  for  fixed  Ng  (propagation  move).  In  order  to  ensure  that  the  Markov 
chain  is  reversible,  the  Gibbs  and  M-H  samplers  of  Equations  (31)  and  (33)  are  used  sys¬ 
tematically  to  update  each  parameter  in  0^  and  0“^,  respectively,  in  the  forward  direction 
for  even  iterations  and  in  the  reverse  direction  for  odd  iterations.  In  other  words,  for  odd 
iterations,  the  procedure  for  updating  0^  in  Equation  (31)  is  applied  backwards  in  the  order 
k  =  Ng,  Ng  —  1, . . .  ,1  (and,  similarly,  for  the  update  of  9^  for  odd  iterations). 

The  propagation  move  described  here  only  involved  the  update  of  0 for  a  fixed  number  of 
sources  Ng.  In  order  to  move  between  configurations  involving  different  numbers  of  sources, 
we  need  to  use  a  reversible-jump  MGMG  procedure  that  allows  updates  between  states  of 
different  dimensions  in  the  hypothesis  space  (e.g.,  such  as  those  associated  with  the  creation 
of  a  new  source  or  the  annihilation  of  an  existing  source) . 

Trans-dimensional  jumps  and  reversible-jump  MCMC 

A  reversible-jump  MGMG  sampling  algorithm  (which  we  use  to  accomodate  between-model 
moves  such  as,  for  example,  a  change  in  the  number  of  sources)  was  first  introduced  by 
Green  [24]  in  the  context  of  the  development  of  a  methodology  for  addressing  the  model 
selection  problem.  The  reversible-jump  MGMG  algorithm  is  very  appealing  in  that  it 
can  be  considered  to  be  a  natural  generalization  of  the  standard  MGMG  algorithm,  with 
the  generalization  allowing  not  only  moves  in  a  parameter  space  of  fixed  dimension  but 
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also  “jumps”  between  model  spaces  Mjv^  of  different  dimensions  (viz.,  involving  different 
numbers  of  sources  Ng)-  More  specifically,  suppose  a  dimension-changing  (or,  between- 
model)  move  of  type  m  is  proposed,  and  the  new  state  0'  is  generated  by  a  deterministic 
invertible  function  g{Q,  v),  where  v  is  a  random  vector  with  distribution  /(v)  [so,  v  /(v)]. 
In  other  words,  g  is  an  operator  that  maps  state  0  and  the  random  vector  v  into  the  new 
state  0h  Green  [24]  demonstrated  that  the  acceptance  probability  of  this  proposed  between- 
model  move  m  has  the  following  form: 


a(0,  0')  =  min 


’  p(0|D,/)r^(0)/(v) 


dg{&,^)  1 

5(0, v)  i’ 


(35) 


where  m'  denotes  the  reverse  move  to  m  and  rm{&)  denotes  the  probability  of  choosing  a 
move  of  type  m  in  the  state  0.  The  final  term  in  the  ratio  of  Equation  (35)  is  the  Jacobian 
which  results  from  the  change  in  variables  associated  with  the  dimension-changing  move. 


For  our  current  problem,  we  consider  two  types  of  dimension-changing  moves:  namely,  a 
creation  move  that  results  in  the  addition  of  a  single  new  source  to  the  current  source 
distribution,  and  an  annihilation  move  that  results  in  the  removal  of  a  single  existing  source 
from  the  current  source  distribution.  To  be  more  specific,  a  creation  operator  C  generates  a 
between-model  move  from  {Ng,  9nJ  ^  in  model  to  (A^^+l,  ^Afa-i-i)  £ 

in  model  so  0'  =  {Ng  +  l,6]\[^+i)  =  C(0)  =  C[{Ng,6]\[^)y  Similarly,  the  reverse 

move,  associated  with  the  annihilation  operator  ,  results  in  an  existing  source  being 
deleted:  so,  0'  =  {Ng  -  1,9n,-i)  =  ct(0)  =  ct  ((iV^,  ^atJ)  where  {Ng  -  1,9n,-i)  G 

^;^l+6A^s 


To  move  from  model  Mjy,  to  using  the  creation  operator  C,  we  propose  the  generation 

of  a  new  source  at  a  location  with  source  strength  Qat^+i  ,  emitting  material  between 

the  activation  and  deactivation  times  and  respectively.  The  “coordinates” 

of  the  new  source  are  generated  by  drawing  random  samples  from  some  proposal  density 
which  we  choose  to  be  the  prior  density  for  each  coordinate:  so,  from  Equations  (19)  to 
(22)  we  sample  yig,N,+i  ~  jl),  Qn,+i  ~  p{Qn,+i\I),  ~  p{T^‘’^^\I),  and 

T^s+i  ^  p{T^‘~^^\T^‘'~^^ ,  I).  Now  assemble  the  parameters  describing  the  new  source  as 
i’  =  {-^g,N,+i,T^‘’~^^,T^‘~'~^,QN,+i),  so  ^AT^+i  =  C{9n,)  =  {9n,,'iP)  for  random  vector 
■0  G  whose  components  are  sampled  as  described  above. 


Suppose  that  at  any  iteration,  we  propose  a  creation  move  with  probability  pc,  with  the 
reverse  (annihilation)  move  having  the  probability  p^^.  The  proposal  probability  for  the 
creation  move  is  then 


=  PeX 


9c(e,0O  =pc  xp(x,,^,+i|/)p(r,^»+V)Mri"“+'|r^+\/)p(Qiv,+i|/) 

1 


^Ns+l^j^Ns+l 

1 


E(P)T^ax(r_.-r,^»+i) 


(36) 


The  probability  (7C|(0^0)  for  the  reverse  (annihilation)  move  is  equal  to  the  probability 
of  choosing  this  move  {p^^)  times  the  probability  of  picking  a  particular  source  from  the 
Ng  +  1  available  sources  for  annihilation,  so 


cT  iV,  +  1 


(37) 
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In  view  of  Equation  (25),  the  ratio  of  posterior  distributions  in  Equation  (35)  for  0'  =  C(0) 
has  the  following  form: 


p(0|D,/) 


=  exp  (^-^(x^(0O  -X^(0))) 

^^p{N,  +  l\I)  11  1  1 

p{N,\I)  E(P)T^ax(r^ax-r,^“  +  ')Qm 


(38) 


Einally,  if  we  substitute  Equations  (36),  (37),  and  (38)  in  Equation  (35),  and  note  that  for 
the  creation  operator  C  the  Jacobian  term  in  Equation  (35)  is  simply  unity,  the  acceptance 
probability  a(0,  0')  for  the  creation  move  0'  =  C(0)  is  given  by 


a(0,  0')  =  min 


p{N,  +  l\I)Pc]  1 
p{N.s\I)  Pc  N.s  +  1 


exp  (x^(0')  -  X^(0)))  I  • 


(39) 


As  demonstrated  by  Green  [24],  a  sufficient  condition  to  ensure  reversibility  of  the  trans- 
dimensional  Markov  chain  is  for  the  acceptance  ratio  for  the  reverse  move  to  be  given  by  the 
reciprocal  of  that  for  the  forward  move.  In  the  present  context,  the  condition  of  detailed 
balance  required  for  Markov  chain  reversibility  implies  that  the  acceptance  probability, 
associated  with  the  annihilation  move  0'  =  ct(0)  for  the  removal  of  a  source  (from  a 
source  distribution  containing  Ng  sources),  must  have  the  following  form: 


q;(0,  0')  =  min 


p{Ng  -1\I)  PC 
p{Ng\I) 


(-i(y{e')-,y(e)))}. 


(40) 


It  only  remains  now  to  specify  the  probabilities  pc  and  p  ^  for  creation  and  annihilation 
moves,  respectively.  0f  course,  with  probabilities  specfied  for  the  creation  and  annihilation 
moves,  the  probability  p-p  for  the  remaining  propagation  move  (source  parameters  updated 
for  a  hypothesis  space  of  fixed  dimension)  is  determined  as  pp  =  1  —pc  ~P^^-  Furthermore, 
it  is  useful  to  allow  the  probabilities  pc  and  to  depend  on  the  number  of  sources  Ng 


in  the  current  state.  To  indicate  explicitly  this  dependence  on  Ng  we  will  augment  the 
notation,  and  express  the  probability  of  a  creation  and  annihilation  move  as  and  p^l , 

Cl 

respectively.  For  Ng  =  1,  p^f  =  0  because  at  least  one  source  must  be  responsible  for 

cT 

the  concentration  measured  by  the  array  of  sensors.  Also,  for  Ng  =  A^s,max)  Pc"  —  0  for 
otherwise  the  preassigned  maximum  number  of  sources  that  may  be  responsible  for  the 
measured  concentration  will  be  exceeded.  For  all  other  cases,  the  probabilities  for  creation 
and  annihilation  moves  will  be  specified  as  follows: 


P 


P^c" 


Ns  +  l 


=  —  mm 
2 

1 

=  —  min  <  1 


p{Ng  +  l\I)\ 


p{Ng\I) 
p{Ng\I) 
^’p{Ng  +  l\I) 


7’ 


(41a) 

(416) 


This  particular  choice  for  and  p  guarantees  that  p^^p{Ns\I)  =  p  p{Ns  +  1|/) 

Cl  cT 

and  ensures  that  the  probability  of  a  dimension-changing  move  at  each  iteration  is  between 
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0.5  and  1  (Green  [24]).  Finally,  it  should  be  noted  that  with  the  dependence  of  pc 
on  Ng  as  in  Equation  (41),  the  ratios  of  these  two  probabilities  in  Equations  (39) 

(40)  need  to  be  interpreted  as  follows:  namely,  p  f/pc  / Pc"  ™  Equation  (39) 

^  '  C  I 

Pc/p^f^Pc‘ 


in  Equation  (40). 


and 

and 

and 


Parallel  tempering  and  Metropolis-coupled  MCMC 

To  overcome  problems  associated  with  a  Markov  chain  that  is  slowly  mixing  in  the  hy¬ 
pothesis  space  and  which  can  cause  the  chain  to  be  stuck  in  a  local  minimum,  we  have 
implemented  a  form  of  parallel  tempering  to  obtain  a  Metropolis-coupled  MCMC  algorithm 
that  improves  the  speed  at  which  the  hypothesis  space  is  explored.  The  Metropolis-coupled 
MCMC  algorithm  has  been  described  by  Geyer  [25]  and  is  based  on  the  idea  of  using  a 
series  of  transition  densities  Ti(0'j0),  T2(0'j0), . .  .,Tj.(0'j0)  with  corresponding  (unnor¬ 
malized)  invariant  (stationary)  distributions  pi(0jD, /),p2(0|D,  .f),  •  •  •  )Pr(0|D,  .f)-  In  our 
current  application  of  Metropolis-coupled  MCMC,  we  run  r  Markov  chains  in  parallel  to 
give  samples  0j*j^  (i  =  1, 2, . . . ,  r),  where  chain  i  has  the  stationary  distribution  given  by 

p,(0lD,/)=p(0lD,/)/«(Dl0,/),  z  =  l,2,...,r,  (42) 

where  G  [0, 1]  (i  =  1, 2, . . . ,  r)  is  an  increasing  sequence  (viz.,  <  Xj  for  i  <  j)  with 
Ai  =  0  and  A,.  =  1. 


The  parameter  A^  in  Equation  (42)  is  interpreted  as  a  tempering  parameter.  This  param¬ 
eter  is  used  to  raise  the  likelihood  function  to  the  A^  power  to  give  a  modified  posterior 
proportional  to  p(0j/)p^*(Dj0, 1).  It  is  noted  that  A^  =  1  is  associated  with  the  desired 
target  posterior  probability  distribution  [cf.  Equation  (25)]  that  we  want  to  sample  from. 
The  other  simulations  correspond  to  a  “ladder”  of  modified  posterior  distributions  indexed 
by  i.  As  the  tempering  parameter  increases  from  zero  to  one,  the  effects  of  the  concen¬ 
tration  data  are  introduced  slowly  through  the  “softened”  likelihood  function  p^*(Dj0,/) 
[Ai  G  [0,1)  for  i  =  1, 2, . . . ,  r  —  1].  In  particular,  A^  <  1  implies  that  the  corresponding 
modified  posterior  distribution  is  broader  (or  flatter)  than  the  actual  posterior  distribution, 
allowing  the  states  to  move  more  freely  in  the  hypothesis  space. 


Let  us  denote  the  state  of  the  i-th  Markov  chain  at  iteration  t  by  0j*j^  (i  =  1, 2, . . . ,  r).  In 
this  form  of  Metropolis-coupled  MCMC  using  parallel  tempering,  the  r  Markov  chains  are 
run  simultaneously  and  at  each  iteration  there  is  a  prescribed  probability  (or,  equivalently, 
on  average  every  A"iter  iterations)  that  a  pair  of  adjacent  chains  (say,  i  and  i  -|-  1  with 
i  =  l,2,...,r  —  l)is  randomly  selected,  and  a  proposal  is  made  to  swap  the  states  of  these 
two  chains.  More  specifically,  if  at  iteration  t,  a  proposed  swapping  operation  between 
chains  i  and  i  -|-  1  is  accepted,  then  we  swap  the  states  of  the  chains  0j*j^  ^ 

^  with  an  acceptance  probability  for  this  swap  given  by 


aswap  =  min 


^  p,(0[^,]lD,/)p,+i(0[yiD,/)| 

’p,(0[JlD,/)p,+i(0«,,lD,/)J 


(43) 


This  swap  enables  the  exchange  of  information  across  the  population  of  r  parallel  simula¬ 
tions.  More  specifically,  this  allows  the  chain  associated  with  the  desired  target  posterior 
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distribution  (viz.,  that  corresponding  to  A,.  =  1)  to  sample  from  remote  regions  of  the 
posterior  distribution,  which  in  turns  facilitates  chain  mobility  in  the  hypothesis  space  and 
ensures  a  more  reliable  relaxation  of  the  chain  into  the  (potentially  exponentially  tiny) 
regions  of  this  space  where  the  posterior  probability  mass  is  concentrated.  We  apply  the 
Metropolis-coupled  MCMC  algorithm  described  here  with  r  =  21.  The  tempering  parame¬ 
ters  Xi  used  on  this  ladder  of  parallel  simulations  is  uniformly  spaced  between  0  and  1. 


Tests  with  synthetic  concentration  data 


In  this  section,  we  present  the  results  of  the  proposed  algorithm  for  multiple  source  re¬ 
construction  for  two  cases:  namely,  case  1  involves  two  unknown  sources,  whereas  case  2 
involves  three  unknown  sourcess.  In  each  of  these  cases,  simulated  concentration  data  will 
be  generated  and  the  Bayesian  inference  scheme  described  above  will  be  applied  to  recon¬ 
struct  the  unknown  source  parameters.  It  should  be  stressed  that  currently  there  are  no 
available  real  concentration  data  sets  that  can  be  used  to  test  our  proposed  algorithm  for  an 
unknown  number  of  multiple  sources;  hence,  the  necessity  to  use  simulated  concentration 
data  sets  instead.  However,  it  should  be  noted  that  to  address  this  deficiency.  Technical 
Panel  9  (TP-9)  of  The  Technical  Cooperation  Program  (TTCP)  Chemical,  Biological  and 
Radiological  Defense  (CBD)  Group  has  proposed  the  design  and  implementation  of  a  coop¬ 
erative  FUsing  Sensor  Information  from  Observing  Networks  (FUSION)  Field  Trial  2007 
(FFT  07).  This  field  experiment  will  be  conducted  at  Tower  Grid  on  U.S.  Army  Dugway 
Proving  Ground  (DPG)  in  September  2007.  The  objective  of  this  field  experiment  will  be 
to  provide  a  comprehensive  tracer  dispersion  and  meteorological  dataset  suitable  for  testing 
current  and  future  chemical  and  biological  sensor  data  fusion  algorithms. 


In  view  of  this,  we  will  simulate  concentration  data  sets  that  conform  with  the  proposed 
design  of  FFT  07.  It  should  be  stressed  that  the  terrain  at  Tower  Grid  is  flat  and  smooth 
with  a  uniform  upwind  fetch  of  more  than  20  km.  The  vegetation  at  this  site  consists 
of  either  widely  spaced  low-lying  shrub  (less  than  0.5  m  in  height)  interspersed  with  bare 
ground,  or  very  short  grass.  The  aerodynamic  roughness  length,  zq,  for  the  site  is  about 
0.015  m  (Yee  et  al.  [26]).  Owing  to  the  horizontal  homogeneity  of  the  site,  the  mean 
wind  flow  and  turbulence  statistics  will  be  assumed  to  be  horizontally  homogeneous  and 
stationary.  Given  the  short  times  and  distances  over  which  we  will  be  modeling  dispersion, 
this  assumption  is  quite  acceptable  (and,  indeed,  at  these  distances  the  dispersion  is  assumed 
to  occur  entirely  within  the  atmospheric  surface  layer  only) .  To  simulate  the  concentration 
seen  by  the  array  of  sensors,  we  use  the  forward-time  LS  model  given  by  Equations  (3) 
and  (5)  with  the  choice  of  the  Kolmogorov  constant  Cq  =  4.8,  the  latter  value  having 
been  obtained  by  calibrating  the  LS  model  against  concentration  data  measured  during 
a  benchmark  field  experiment  Project  Prairie  Grass  (Wilson  et  al.  [27]).  The  wind  flow 
statistics  required  by  the  LS  model  are  prescribed  as  follows  in  accordance  primarily  with 
well-known  surface-layer  relations  based  on  Monin-Obukhov  theory  (Stull  [28]). 


Let  X,  y,  and  z  represent  coordinates  in  the  alongwind,  crosswind,  and  vertical  directions, 
respectively.  The  instantaneous  (mean)  velocities  in  these  three  directions  will  be  denoted 
by  u  ([/),  V  (U),  and  w  (IT),  respectively.  Furthermore,  u' ,  v' ,  and  w'  are  departures  of  u, 
V,  and  w  from  its  mean  values  U,  V,  and  W,  respectively.  Then,  the  stability-dependent 
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mean  wind  speed  profile  in  the  atmospheric  surface  layer  has  the  following  form: 


U{z)  = 


(44) 


where  is  the  friction  velocity,  k-u  ~  0.4  is  von  Karman’s  constant,  L  is  the  Obukhov 
length  scale  and  V’  is  the  stability  function.  For  unstable  stratification  (L  <  0),  V’  assumes 
the  form 


2arctan(x)  + 


where 


and,  for  stable  stratification  (L  >  0) 


L  <  0, 


Note  that  the  mean  wind  direction  is  assumed  to  be  along  the  x-axis,  so  14  =  0. 


(45a) 


(456) 


(45c) 


The  parametrizations  we  use  for  the  turbulence  velocity  flow  statistics  follow  those  suggested 
by  Rodean  [29] .  For  the  coordinate  system  used  here  where  the  x-axis  is  chosen  to  lie  along 
the  mean  wind  direction,  the  velocity  covariances  u'v'  and  v'w'  are  zero  in  the  atmospheric 
surface  layer;  and,  the  only  remaining  velocity  covariance  u'w'  is  assumed  to  be  constant 
and  equal  to  —u1  (by  definition).  The  alongwind  and  crosswind  velocity  variances  and 
cr^,  respectively,  are  given  by 
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where  H  is  the  boundary-layer  height.  The  vertical  velocity  variance  <7^  is  given  by 
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Finally,  the  dissipation  e  of  turbulence  kinetic  energy  is  parameterized  as  follows: 
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For  the  two  examples  used  to  test  the  Bayesian  inference  algorithm  for  source  reconstruc¬ 
tion,  the  concentration  measured  by  a  particular  sensor  in  the  array  will  result  from  the 
superposition  of  two  or  more  partially  overlapping  plumes  produced  by  continuously  emit¬ 
ting  sources.  It  is  noted  that  for  multiple  continuous  point  sources,  ^  — oo  and  ^  oo 
[A:  =  1,2  . .  .,Ns]  cf.  Equations  (13)  and  (14)].  In  consequence,  in  this  case,  the  only  rel¬ 
evant  parameters  are  the  source  location  ^  and  the  emission  rate  Qk  for  each  source 
{k  =  1,  2, . . . ,  A^s).  Furthermore,  it  is  noted  that  the  concentration  field  C  and  adjunct 
concentration  field  C*  are  independent  of  time  [viz.,  (^(x,  t)  =  C'(x)  and  t)  =  (^^(x)]. 

Finally,  for  the  simulations,  it  is  assumed  that  all  the  sources  are  emitting  at  ground  level 
{z  =  0)  and  that  this  is  known  a  priori  (viz.,  this  knowledge  is  considered  to  be  part  of  the 
background  information  I).  As  a  result,  the  unknown  location  parameters  for  each  source 
are  its  alongwind  (xg)  and  crosswind  (y^)  positions,  only. 


Example  1:  two  unknown  sources 


For  the  first  example,  we  synthesized  artificial  concentration  data  for  the  case  of  two  sources 
which  were  located  upwind  of  an  array  consisting  of  42  detectors  arranged  as  shown  in 
Figure  1.  Each  detector  was  placed  at  a  height  of  1.5  m  above  ground  level.  The  two 
ground  level  continuously  emitting  sources  were  located  at  {xg^yg)  =  (—50.0,0.0)  m  and 
at  (—250.0,0.0)  m  and  each  source  had  an  emission  rate  oi  Q  =  Qg  =  1.0  g  s“^.  The 
flow  statistics  for  the  atmospheric  surface  layer  (through  which  the  dispersion  occurred) 
require  three  surface  layer  parameters:  tt*,  L,  and  zq  which  for  the  simulations  assumed  the 
following  values:  =  0.25  m  s“^,  L  =  10000  m,  and  zq  =  0.015  m  (roughness  height  for 

Tower  Grid  at  U.S.  Army  Dugway  Proving  Ground).  It  is  assumed  that  the  boundary-layer 
height  was  H  =  1000  m.  The  concentration  data  generated  using  the  forward-time  LS  model 
(with  these  flow  statistics  as  input)  were  embedded  within  white  and  normally  distributed 
noise  with  a  standard  deviation  equal  to  10%  of  the  true  concentration  amplitude. 


We  applied  our  proposed  algorithm  for  source  reconstruction  to  the  simulated  concentra¬ 
tion  data.  We  perform  Metropolis-coupled  MGMG  sampling  with  r  =  21  chains  with  swaps 
between  these  chains  attempted  at  every  A"iter  =  25  iterations  on  average.  The  proposed 
algorithm  randomly  initializes  all  unknown  parameters  (source  location,  emission  rate)  in 
accordance  to  the  prescribed  prior  distributions  for  these  parameters  [cf.  Equations  (19) 
and  (20)1  with  Qmax  =  100  kg  s“^  and  V  =  [-2000,-25]  m  x  [—500,500]  m  providing 
the  prior  bounds  on  the  emission  rate  Qg  and  on  the  source  location  {xg,yg),  respectively. 
The  initial  number  of  sources  was  randomly  assigned  in  accordance  to  the  prior  distribu¬ 
tion  for  Ng  given  by  Equation  (18),  where  A^s,max  =  4  (maximum  allowable  number  of 
sources)  and  the  hyperparameter  F  (expected  number  of  sources)  is  initialized  to  1  (hence, 
the  prior  distribution  favors  the  wrong  choice  for  the  number  of  sources).  For  this  example, 
the  MGMG  algorithm  was  run  for  100,000  iterations,  with  the  first  50,000  iterations  corre¬ 
sponding  (conservatively)  to  the  burn-in  with  the  result  that  these  samples  were  discarded 
from  the  analysis.  The  remaining  50,000  post  burn-in  samples  were  used  for  the  posterior 
inference. 


Figure  2  (top)  displays  changes  in  Ng  (number  of  sources)  as  a  function  of  the  iteration 
number  for  the  first  500  post  burn-in  iterations  of  the  reversible-jump  MGMG  sampler. 
Note  the  dimension-changing  moves  involving  creation  of  a  new  source  (e.g.,  transitions 
from  Ng=2\,oNg  =  2>OT  from  Ng  =  3  to  =  4)  or  annihilation  of  an  existing  source 
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(e.g.,  transitions  from  =  4  to  =  3  or  from  Ng  =  3  to  Ng  =  2).  Interestingly,  after 
convergence  of  the  Markov  chain  to  the  stationary  distribution,  moves  across  models  from 
A^s  =  2  to  A^s  =  1  (and,  vice-versa)  do  not  occur.  Indeed,  Figure  2  (bottom)  which  exhibits 
the  probability  distribution,  p{Ng)  =  p(A^s|D,  I),  for  the  number  of  sources  estimated  from 
the  50,000  post  burn-in  samples,  shows  that  the  hypothesis  Ng  =  1  source  is  excluded. 
The  simulations  settle  in  a  distribution  which  favors  equally  (approximately  or  better)  the 
hypotheses  Ng  =  2  or  3,  with  a  smaller  probability  for  A^^  =  4. 

Figure  3  exhibits  histograms  of  the  alongwind  location  Xg  (top  panel),  crosswind  location 
i/g  (middle  panel) ,  and  emission  rate  Qg  (bottom  panel) .  These  histograms  were  constructed 
from  the  subset  of  post  burn-in  samples  consisting  of  two  source  objects  (i.e.,  samples  with 
Ng  =  2).  The  bimodality  of  the  histogram  of  Xg  with  mode  locations  at  Xg  ~  —50  m 
and  Xg  ~  —250  m  indicates  that  the  alongwind  positions  of  the  two  sources  have  been 
correctly  identified.  Because  both  sources  are  located  at  =  0  m  and  have  an  emission 
rate  (7^  =  1  g  s“^,  the  single  mode  in  the  histograms  at  ~  0  m  and  (7^  ~  1-0  g  s“^  indicates 
that  these  source  parameters  have  been  correctly  inferred  for  the  two  sources.  The  samples 
corresponding  to  each  mode  in  the  histogram  for  Xg  can  be  extracted  separately.  Figures  4 
and  5  display  the  first  12,000  samples  obtained  from  each  of  these  two  modes  in  the  form 
of  trace  histories  of  the  source  parameters  {xg,yg,qg)  corresponding  to  each  mode.  It  can 
be  seen  that  the  source  parameters  associated  with  each  mode  (and,  corresponding  to  a 
particular  source  object)  are  estimated  well.  Finally,  it  should  be  stressed  that  labels  used 
for  the  source  objects  in  Figures  4  and  5  [e.g.,  source  object  1  is  identified  to  be  the  source 
at  {xg,yg)  =  (—250,0)  m,  whereas  source  object  2  is  the  source  at  {xg,yg)  =  (—50,0)  m] 
are  completely  arbitrary  here  because  the  posterior  probability  distribution  of  the  source 
parameters  [cf.  Equation  (25)]  is  invariant  under  a  reordering  (relabelling)  of  the  identifiers 
used  for  each  source  object.  This  degeneracy  simply  corresponds  to  different  (but  equivalent) 
identifications  of  what  is  meant  by  source  1,  source  2,  etc. 

Figures  6  and  7  display  histograms  of  the  alongwind  location  Xg  (top  panel),  crosswind 
location  y^  (middle  panel),  and  emission  rate  qg  (bottom  panel)  which  were  obtained,  re¬ 
spectively,  from  the  subset  of  post  burn-in  samples  consisting  of  three  {Ng  =  3)  and  four 
{Ng  =  4)  source  objects.  Note  from  both  these  figures  that  there  really  only  exists  two 
modes  in  the  histograms  of  Xg,  and  it  is  seen  that  these  modes  correspond  to  the  correct 
alongwind  locations  of  the  two  sources.  Furthermore,  the  extra  source  object  in  the  his¬ 
tograms  in  Figure  6  and  the  two  extra  source  objects  in  the  histograms  of  Figure  7  are 
randomly  distributed  in  D  (viz.,  are  not  concentrated  into  a  significant  cluster  or  clusters 
of  points).  Interestingly,  these  extra  source  objects  have  emission  rates  associated  with  the 
mode  in  the  histogram  of  at  0  g  s“^.  The  mode  in  the  histogram  of  (7^,  at  1.0  g  s“^ 
is  associated  with  the  modes  in  the  histogram  of  Xg  at  ~  —250  m  and  ~  —50  m  and  the 
mode  in  the  histogram  of  y^  at  ~  0  m.  Hence,  even  the  source  reconstruction  based  on 
samples  with  Ng  =  3  or  Ng  =  4  only  give  two  main  clusters  of  points  in  D  located  at  the 
true  positions  of  the  two  sources  and  these  two  main  clusters  are  associated  with  the  mode 
in  the  histogram  of  qg  at  ~  1.0  g  s“^  (corresponding  to  the  true  emission  rate  of  the  two 
sources). 

From  Figures  3,  6,  and  7  corresponding  to  the  case  Ng  =  2,  3  and  4,  respectively,  we  see 
that  there  exists  two  clusters  in  which  the  samples  are  concentrated.  These  two  clusters 
are  associated  with  two  source  objects.  We  have  obtained  all  samples  for  Ng  =  2,  3,  and  4 
for  each  of  these  two  clusters  and  plotted  histograms  of  the  parameters  y^,  (7s}  for  each 
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detected  source  object  in  Figures  8  and  9.  The  mean,  standard  deviation,  and  lower  and 
upper  bounds  for  the  95%  highest  posterior  density  (HPD)  interval  of  the  parameters  for 
each  detected  source  object  were  calculated  from  the  samples  for  Ns  =  2,  3,  and  4  contained 
in  each  of  the  two  identified  clusters  and  the  results  are  summarized  in  Table  1.  From  this 
information,  we  see  that  the  two  source  objects  have  been  characterized  to  a  reasonable 
accuracy. 

Example  2:  three  unknown  sources 

For  our  second  example,  we  synthesized  artificial  concentration  data  using  a  forward-time 
LS  model  for  the  case  of  three  sources  which  were  located  upwind  of  an  array  consisting  of 
66  detectors  arranged  as  shown  in  Figure  10.  Each  detector  was  placed  at  a  height  of  1.5  m 
above  ground  level.  The  three  ground  level  continuously  emitting  sources  were  located  at 
{xs,ys)  =  (—50.0,0.0)  m,  (-125.0,-20.0)  m  and  (—300.0,20.0)  m  and  each  source  had 
an  emission  rate  of  Q  =  (7^  =  1.0  g  s“^.  The  properties  of  the  atmospheric  surface  layer, 
through  which  the  dispersion  from  the  three  sources  occurred,  was  the  same  as  that  for  the 
first  example  above.  The  synthetic  concentration  data  generated  using  the  forward-time  LS 
model  were  corrupted  with  white  and  normally  distributed  noise  with  a  standard  deviation 
equal  to  10%  of  the  true  concentration  amplitude. 

We  applied  our  proposed  algorithm  for  source  reconstruction  to  the  simulated  concentration 
data.  The  Metropolis-coupled  reversible-jump  MCMC  algorithm  used  the  same  parameters 
as  described  previously  for  Example  1,  with  the  following  exceptions.  Rather  than  specifying 
the  prior  distribution  for  Q  as  in  Equation  (19),  we  have  used  instead 

p(Qfc|L)  =  (1 -7)'5(Qfc) +  7^(0,Q,„ax)(Qfc)/Qmax,  k  =  1 ,2 , . . . ,  Ns,  (50) 

where  7  is  an  intermittency  factor,  defined  as  the  probability  that  the  source  is  turned  on 
(viz.,  Qk  >  0).  Equation  (50),  unlike  Equation  (19),  allows  for  the  hypothesis  that  any 
given  source  in  the  domain  can  be  turned  off  {Qk  =  0)  with  a  finite  probability  (=1  —  7). 
Eor  the  source  reconstruction  here,  we  have  chosen  the  hyperparameter  7  =  0.2.  The  initial 
number  of  sources  was  randomly  assigned  in  accordance  to  the  prior  distribution  for  Ns 
given  by  Equation  (18),  where  W,max  =  8  (maximum  allowable  number  of  sources)  and 
the  hyperparameter  F  (expected  number  of  sources)  is  initialized  to  1.  In  consequence, 
our  initial  specification  for  the  prior  distribution  for  Ns  favors  the  wrong  choice  for  the 
actual  number  of  sources.  We  note  that  the  MCMC  sampler  moved  quickly  from  the  low 
likelihood  of  the  starting  point  to  a  region  in  the  hypothesis  space  with  higher  likelihood. 
As  in  Example  1,  the  MCMC  algorithm  was  run  for  100,000  iterations,  with  the  first  50,000 
iterations  corresponding  to  the  burn-in  with  the  result  that  these  samples  were  discarded 
from  the  analysis.  The  remaining  50,000  post  burn-in  samples  were  used  for  the  posterior 
inference,  the  results  of  which  will  be  presented  below. 

An  estimate  of  the  posterior  distribution  for  Ns,  obtained  from  a  histogram  of  the  number  of 
samples  obtained  in  each  subspace  of  different  dimension,  is  exhibited  in  Eigure  11.  We  note 
that  the  most  probable  number  of  source  objects  is  Ns  =  3,  which  corresponds  to  the  correct 
number  of  sources.  This  result  is  obtained  in  spite  of  the  fact  that  the  prior  distribution 
for  Ns  was  initialized  with  F  =  1  (a  priori  expected  number  of  sources),  implying  that 
the  algorithm  is  not  sensitive  to  this  hyperparameter.  The  information  embodied  in  the 
concentration  data  was  sufficient  to  move  the  simulations  towards  the  more  complex  model 
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with  Ng  =  3  sources.  Interestingly,  it  was  found  that  with  the  prior  specification  for  Qk  given 
by  Equation  (50),  many  samples  for  more  complex  models  involving  Ng  >  3  sources  had  the 
additional  sources  turned  off  (viz.,  =0).  It  should  be  noted  that  for  samples  involving 
Ng  sources,  but  with  of  these  sources  turned  off,  the  source  model  here  is  considered 
to  have  Ng  —  Nd  sources,  and  not  Ng  (owing  to  the  fact  that  sources  that  are  turned  off 
do  not  contribute  to  the  model  concentration  seen  at  the  detectors).  In  consequence,  some 
of  the  samples  corresponding  to  the  mode  of  p{Ng)  =  p{Ng\'D,I)  in  Figure  11  at  Ng  =  3 
corresponds  to  more  complex  models  involving  more  than  3  sources,  but  with  the  additional 
sources  turned  off. 

Figure  12  displays  histograms  of  the  source  parameters  {xg,yg,qg}  constructed  from  all 
samples  corresponding  to  Ng  =  3  (mode  for  p{Ng)  in  Figure  11).  Furthermore,  we  plot  the 
samples  obtained  for  the  case  Ng  =  3,  projected  onto  the  (x^,  ys)-subspace  in  Figure  13. 
Both  Figures  12  and  13  indicate  that  samples  tend  to  cluster  in  three  main  regions  of 
the  hypothesis  space,  which  can  be  identified  with  three  source  objects.  Histograms  of 
the  parameters  associated  with  each  of  these  three  source  objects  are  exhibited  in  Figures 
14,  15  and  16.  The  mean,  standard  deviation,  and  lower  and  upper  bounds  for  the  95% 
HPD  interval  of  the  source  parameters  for  each  identified  source  object  are  summarized  in 
Table  2.  Comparing  these  values  of  the  parameters  of  the  recovered  source  objects  with  the 
parameters  for  the  three  actual  sources,  we  see  that  the  algorithm  has  adequately  recovered 
the  true  parameters  for  each  source  object,  within  the  stated  errors.  Not  surprisingly,  the 
parameters  for  the  source  furtherest  away  from  the  array  of  sensors  were  determined  with 
the  largest  uncertainty  (as  measured  either  by  the  standard  deviation,  or  the  95%  HPD 
interval) . 


Conclusions 


In  this  paper,  we  have  presented  a  Bayesian  inference  approach  for  multiple  source  recon¬ 
struction  from  a  limited  number  of  noisy  concentration  data  obtained  from  an  array  of 
sensors  for  the  case  where  the  number  of  sources  is  unknown  a  priori.  In  this  approach,  we 
use  a  model  that  relates  the  source  distribution  to  the  concentration  data  (source-receptor 
relationship) ,  and  then  apply  Bayesian  probability  theory  to  formulate  the  posterior  density 
function  for  the  source  parameters  including  the  number  of  sources.  The  evaluation  of  the 
posterior  density  function  and  of  its  features  of  interest  requires  a  numerical  procedure.  To 
this  end,  the  computational  algorithm  required  here  is  implemented  using  a  reversible-jump 
Markov  chain  Monte  Carlo  method  to  draw  samples  from  the  posterior  density  function. 
We  showed  how  to  design  creation  and  annihilation  moves  that  allow  the  Markov  chain  to 
jump  between  hypothesis  spaces  corresponding  to  different  numbers  of  sources  in  the  source 
distribution. 

The  new  methodology  has  been  sucessfully  applied  to  simulated  concentration  data  corre¬ 
sponding  to  two  and  three  source  cases.  It  is  shown  that  the  proposed  method  performs 
well:  the  number  of  sources  is  correctly  identified  using  the  procedure,  and  the  parameters 
(e.g.,  location,  emission  rate)  that  characterize  each  identified  source  are  reliably  estimated. 
In  addition,  the  methodology  provides  a  rigorous  determination  of  the  uncertainty  (e.g., 
standard  deviation,  credible  intervals)  in  the  inference  of  the  source  parameters,  hence 
extending  the  potential  of  the  methodology  as  a  tool  for  quantitative  source  reconstruction. 
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TABLE  1 

The  mean,  standard  deviation  (S.D.),  and  lower  and  upper  bounds  of  the  95%  highest 
posterior  density  (HPD)  interval  of  the  parameters  Vs.k  and  Qg^k  {k  =  1,2)  calculated 
from  the  samples  for  each  cluster  in  Figures  3,  6  and  7  for  Ng  =  2,  3  and  4,  respectively. 
These  are  the  parameters  that  characterize  the  two  source  objects  {k  =  1,2)  associated  with 
the  two  clusters  of  samples. 


Parameter 

Mean 

S.D. 

95%  HPD 

k  =  l 

Xg  (m) 

-254 

11 

(-276,-234) 

Vs  (m) 

-0.13 

0.37 

(-0.85,0.58) 

Q.S  (g  s“i) 

1.04 

0.05 

(0.95,1.14) 

k  =  2 

Xg  (m) 

-50 

0.2 

(-50.5,-49.6) 

Vs  (m) 

-0.00195 

0.036 

(-0.07,0.07) 

q.s  (g  S“^) 

1.02 

0.01 

(0.99,1.05) 
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TABLE  2 

The  mean,  standard  deviation  (S.D.),  and  lower  and  upper  bounds  of  the  95%  highest 
posterior  density  (HPD)  interval  of  the  parameters  Vs.k  and  Qg^k  {k  =  1,2,?,)  calculated 
from  the  samples  for  each  cluster  in  Figures  12  and  13  for  Ng  =  3.  These  are  the  parameters 
that  characterize  the  three  source  objects  {k  =  1,2,3)  associated  with  the  three  clusters  of 
samples. 


Parameter 

Mean 

S.D. 

95%  HPD 

k  =  l 

Xg  (m) 

-288 

16 

(-323,-260) 

Vs  (m) 

20.8 

1.1 

(18.4,22.6) 

<ls  (g  s“^) 

0.956 

0.094 

(0.80,1.16) 

k  =  2 

Xg  (m) 

-124 

1.7 

(-127.5,-121.3) 

Vs  (m) 

-20.1 

0.2 

(-20.5,  -19.7) 

<ls  (g  s“^) 

1.03 

0.03 

(0.97,1.09) 

k  =  3 

Xg  (m) 

-50.1 

0.2 

(-50.5,-49.6) 

Vs  (m) 

-0.0073 

0.039 

(-0.084,0.067) 

<ls  (g  s“^) 

1.03 

0.02 

(0.99,1.05) 
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Source-Detector  Configuration 


Figure  1,  Two  continuous  point  sources  located  upwind  of  an  array  of  42  deteetors  arranged 
as  shown.  A  solid  dot  denotes  a  source  and  a  solid  square  denotes  a  sensor. 


DRDC  Suffield  TR  2007-085 


29 


Instantaneous  estimates  of  N  (first  500  iterations) 
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Figure  2,  Instantaneous  estimates  of  Ng  for  the  first  500  post  burn-in  iterations  of  the 
reversible  jump  MCMC  sampler  (top)  and  the  posterior  probability  distribution  for  the 
number  of  sources  (bottom)  estimated  using  the  50,000  post  bum-in  samples. 
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Histogram  of  Source  Parameters  (2  Sources) 
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Figure  3,  Histograms  of  the  alongwind  location  (top  panel),  crosswind  location  (middle 
panel),  and  emission  rate  (bottom  panel)  of  the  two  sources  obtained  from  post  bum-in 
samples  with  Ng  =  2.  The  vertical  line(s)  in  each  panel  mark  the  tme  values  for  the  source 
parameters. 
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Figure  4.  Example  of  a  trace  of  source  parameters  for  a  subset  of  12,000  samples  extracted 
from  the  samples  associated  with  the  mode  at  =  -250  m  in  the  histogram  of  exhibited  in 
Figure  3  (top  panel).  The  solid  horizontal  lines  are  the  true  values  for  the  source  parameters. 
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Sourcc-Dctcctor  Configuration 
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Figure  5.  Example  of  a  trace  of  source  parameters  for  a  subset  of  12,000  samples  extracted 
from  the  samples  associated  with  the  mode  at  =  -50  m  in  the  histogram  of  x^  exhibited  in 
Figure  3  (top  panel).  The  solid  horizontal  lines  are  the  true  values  for  the  source  parameters. 
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^  ^o4  Histogram  of  Source  Parameters  (3  Sources) 
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Figure  6,  Histograms  of  the  alongwind  location  (top  panel),  crosswind  location  (middle 
panel),  and  emission  rate  (bottom  panel)  of  the  two  sources  obtained  from  post  bum-in 
samples  with  Ng  =  3.  The  vertical  line(s)  in  each  panel  mark  the  tme  values  for  the  source 
parameters. 
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Figure  7,  Histograms  of  the  alongwind  location  (top  panel),  crosswind  location  (middle 
panel),  and  emission  rate  (bottom  panel)  of  the  two  sources  obtained  from  post  bum-in 
samples  with  Ng  =  4.  The  vertical  line(s)  in  each  panel  mark  the  tme  values  for  the  source 
parameters. 


DRDC  Suffield  TR  2007-085 


35 


Source  1 


y,(m) 


Sourcc-Dctcctor  Configuration 


-400  -300  -200  -100 

r  (m) 


200 


Figure  8,  Histograms  for  the  three  parameters,  namely  alongwind  location  (top  left  frame), 
crosswind  location  (top  right  frame)  and  emission  rate  (bottom  left  frame)  that  characterize 
source  1  (bottom  right  frame).  The  solid  vertical  line  indicates  the  true  value  of  the  parameter 
and  the  dashed  vertical  line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the 
posterior  mean  of  the  associated  marginal  posterior  distribution. 
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Figure  9,  Histograms  for  the  three  parameters,  namely  alongwind  location  (top  left  frame), 
crosswind  location  (top  right  frame)  and  emission  rate  (bottom  left  frame)  that  characterize 
source  2  (bottom  right  frame).  The  solid  vertical  line  indicates  the  true  value  of  the  parameter 
and  the  dashed  vertical  line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the 
posterior  mean  of  the  associated  marginal  posterior  distribution. 
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Source-Detector  Configuration 


Figure  10,  Three  eontinuous  point  sourees  located  upwind  of  an  array  of  66  detectors 
arranged  as  shown.  A  solid  dot  denotes  a  source  and  a  solid  square  denotes  a  sensor. 
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Figure  11,  The  posterior  distribution  for  the  number  of  sources  estimated  using  the  50,000 
post  burn-in  samples. 
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^  ^q4  Histogram  of  Source  Parameters  (3  Sources) 
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Figure  12,  Histograms  of  the  alongwind  location  (top  panel),  crosswind  location  (middle 
panel),  and  emission  rate  (bottom  panel)  of  the  three  sources  obtained  from  post  bum-in 
samples  with  Ng  =  3.  The  vertical  line(s)  in  each  panel  mark  the  tme  values  for  the  source 
parameters. 
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Samples  Drawn  from  Posterior 


Figure  13,  The  post  burn-in  samples  obtained  for  Ng  =  3  projeeted  onto  the  (xg,yg)  subspace. 
The  true  alongwind  and  crosswind  locations  of  the  three  sources  are  indicated  by  the  open 
circles. 
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Figure  14,  Histograms  for  the  three  parameters,  namely  alongwind  location  (top  left  frame), 
crosswind  location  (top  right  frame)  and  emission  rate  (bottom  left  frame)  that  characterize 
source  1  (bottom  right  frame).  The  solid  vertical  line  indicates  the  true  value  of  the  parameter 
and  the  dashed  vertical  line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the 
posterior  mean  of  the  associated  marginal  posterior  distribution. 
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Figure  15,  Histograms  for  the  three  parameters,  namely  alongwind  location  (top  left  frame), 
crosswind  location  (top  right  frame)  and  emission  rate  (bottom  left  frame)  that  characterize 
source  2  (bottom  right  frame).  The  solid  vertical  line  indicates  the  true  value  of  the  parameter 
and  the  dashed  vertical  line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the 
posterior  mean  of  the  associated  marginal  posterior  distribution. 
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Figure  16,  Histograms  for  the  three  parameters,  namely  alongwind  location  (top  left  frame), 
crosswind  location  (top  right  frame)  and  emission  rate  (bottom  left  frame)  that  characterize 
source  3  (bottom  right  frame).  The  solid  vertical  line  indicates  the  true  value  of  the  parameter 
and  the  dashed  vertical  line  corresponds  to  the  best  estimate  of  the  parameter  obtained  as  the 
posterior  mean  of  the  associated  marginal  posterior  distribution. 
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sources  when  the  number  of  sources  is  unknown  a  priori.  The  problem  is  solved  using  a  Bayesian 
probabilistic  inferential  framework  in  which  Bayesian  probability  theory  is  used  to  derive  the  posterior 
probability  density  function  for  the  number  of  sources  and  for  the  parameters  (e.g.,  location,  emission  rate, 
release  duration)  that  characterize  each  source.  A  mapping  (or,  source-receptor  relationship)  that  relates 
a  multiple  source  distribution  to  the  concentration  data  measured  by  the  array  of  sensors  is  formulated 
based  on  a  forward-time  Lagrangian  stochastic  model.  A  computationally  efficient  methodology  for 
determination  of  the  likelihood  function  for  the  problem,  based  on  an  adjoint  representation  of  the  source- 
receptor  relationship  and  realized  in  terms  of  a  backward-time  Lagrangian  stochastic  model,  is  described. 
An  efficient  computational  algorithm  based  on  a  reversible  jump  Markov  chain  Monte  Carlo  (MCMC) 
method  is  formulated  and  implemented  to  draw  samples  from  the  posterior  density  function  of  the  source 
parameters.  This  methodology  allows  the  MCMC  method  to  jump  between  the  hypothesis  spaces 
corresponding  to  different  numbers  of  sources  in  the  source  distribution  and,  thereby,  allows  a  sample 
from  the  full  joint  posterior  distribution  of  the  number  of  sources  and  source  parameters  to  be  obtained. 

The  proposed  methodology  for  source  reconstruction  is  tested  using  synthetic  concentration  data 
generated  for  cases  involving  two  and  three  unknown  sources. 
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